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FOREWORD 


This  report  is  submitted  to  the  Air  Force  Weapons  Laboratory,  Kirtland  Air 
Force  Base,  New  Mexico,  in  accordance  with  the  requirements  set  forth  in  AF 
Contract  AF  29(601)-6061  and  SA  No.  1  to  AF  Contract  AF  29(601)-6014,  as  amen¬ 
ded.  This  report  is  divided  into  two  volumes.  Volume  I  describes  the  theoretical 
development  and  numerical  methods  used  in  the  generation  of  the  computer  codes. 
Volume  II  consists  of  seven  independently  bound  sections,  each  of  which  trea.s  the 
programming,  coding,  and  operation  of  a  specific  computer  code.  They  are  prepared 
as  an  operator’s  manual. 

Results  of  other  work  performed  under  these  contracts,  which  consisted  of 
application  of  existing  codes  to  specific  radiation  analysis,  have  been  previously 
transmitted  to  the  Air  Force  Weapons  Laboratory  and  are  not  included  in  this  report. 

The  work  has  been  done  under  the  joint  supervision  of  J.  A.  Barton  and 
B.  W.  Mar  within  the  Nuclear  and  Space  Physics  Section  under  the  direction  of 
G.  L.  Keister.  The  technical  contributions  were  prepared  by  the  following  authors: 
Section  2,  W.  R.  Doherty;  Section  3,  J.  R.  Benbrook,  W.  R.  Doherty,  W.  R.  Sheldon, 
and  J.  R.  Thomas;  Section  4,  B.  W.  Mar  and  K.  Moriyasu;  Section  5,  B.  W.  Mar; 
Section  6,  M.  C.  Wilkinson;  Section  7,  B.  W.  Mar;  Section  8,  W.  R.  Doherty;  and 
Appendix,  K.  Moriyasu. 


ABSTRACT 


The  objective  of  this  program  is  to  develop  mathematical  representations  and 
computer  codes  for  determining  the  absorbed  dose  and  the  linear-energy-transfer 
spectrum  at  arbitrary  points  within  a  spacecraft  on  an  arbitrary  trajectory  in  space. 
The  radiation  environment/  including  spectral/  angulai;and  temporal  dependence, 
was  represented  mathematically/  and  calculation  of  the  flux  and  flux  rate  of  each 
particle  type  was  coded.  Specifically/  codes  were  devised  to  (1)  determine  the 
spacecraft  trajectory;  (2)  convert  the  trajectory  in  geographic  coordinates  to  mag¬ 
netic  coordinates  B,  L,  and  t  and  to  equivalent  magnetic  dipole  coordinates  R, 

X/  and  t;  (3)  calculate  the  unit-flux  dose  or  total  dose  of  penetrating  and  secondary 
radiation  at  a  dose  point  within  the  spacecraft;  and  (4)  integrate  the  various  subrou¬ 
tines  and  instruct  a  plotter  to  plot  program  output. 

The  computer  program  was  devised  for  maximum  growth  capability  and  flexi¬ 
bility,  but  there  are  inevitable  shortcomings  in  view  of  the  rapid  development  of 
space  research.  Recommendations  are  made  for  further  studies. 
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1.0  INTRODUCTION 


Prior  to  1956/  the  radiation  hazard  to  manned  spaceflight  was  not  apparent. 

The  subsequent  rapid  increase  in  knowledge  of  solar  particle  events  and  the  discovery 
of  geomagnetical ly  trapped  particles  in  1958  has  led  to  the  necessity  of  considering 
space  radiation  in  the  development  of  manned  space  systems.  To  meet  the  need  for 
preliminary  design  information,  digital  computer  programs  were  developed  to  predict 
the  radiation  fluxes  and  radiation  doses  encountered  in  space  missions. 

The  first-generation  programs  computing  space-radiation  doses  were  developed 
to  provide  preliminary  design  information,  but  they  were  based  on  meager  informa¬ 
tion  about  the  radiation  environment  and  on  a  simplified  shielding  analysis.  However, 
these  programs  have  provided  a  basis  for  the  development  of  better  programs,  incor¬ 
porating  new  environmental  data  and  a  more  sophisticated  shielding  analysis. 

The  development  of  the  second-generation  computer  program,  containing  addi¬ 
tions  and  modifications  to  the  original  program  developed  by  The  Boeing  Company, 
has  been  particularly  concerned  with  the  trajectory  and  environment  codes  and  the 
radiation -transport  codes. 

Throughout  the  development  of  this  new  program,  efforts  have  been  made  to 
provide  maximum  flexibility,  growth  potential,  and  minimum  machine  time  compati¬ 
ble  with  detail  and  accuracy.  To  provide  flexibility,  the  computer  codes  are  so 
designed  that  the  analyst  has  the  option  of  using  them  independently.  For  example, 
the  shielding  codes  may  be  used  to  determine  the  optimum  vehicle  structure  and 
arrangement  of  internal  equipment  to  provide  maximum  radiation  shielding;  the  tra¬ 
jectory  code  may  be  used  to  select  an  orbit  or  trajectory  incurring  an  acceptable 
radiation  flux  or  flux  rate;  or  a  variety  of  B-  L  maps  may  be  employed  to  determine 
the  effect  of  the  environmental  model.  To  provide  growth  capability,  the  program 
can  readily  accept  new  experimental  data  on  the  radiation  environment  and  on  the 
interaction  of  radiation  with  matter.  Several  techniques  have  been  employed  to 
reduce  the  computer  machine  time,  but  the  analyst  has  the  option  of  increasing  the 
accuracy  of  the  flux  calculation  at  the  expense  of  increased  machine  time.  With  an 
understanding  of  the  assumptions,  limitations,  optional  modes,  and  ramifications  of 
the  routines,  and  with  the  exercise  of  good  judgment,  the  program  will  provide  valid, 
accurate,  and  economical  solutions  to  space -radiation  problems. 

The  succeeding  sections  include  the  theoretical  development  of  the  mathemati¬ 
cal  representations,  a  discussion  of  the  reliability  and  validity  of  the  assumptions 
that  have  been  made,  and  an  evaluation  of  the  results  obtained  with  these  codes. 

The  sequence  of  these  sections  is  in  parallel  with  the  following  general  methods  of 
determining  the  radiation  fluxes  and  doses  for  space  systems: 
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1)  Calculation  of  the  spacecraft  trajectory  in  B,  L,  and  t  coordinates; 

2)  Devising  a  mathematical  representation  of  the  space -radiation  environment-- 
the  radiation  environment-  includes  geomagnetical ly  trapped  radiation  (Van  Allen 
belts)/  solar  particle  event  radiation/  and  galactic  cosmic  radiation; 

3)  Determination  of  the  radiation  flux  and  energy  spectra  encountered  in  a 
given  space  mission; 

4)  Calculation  of  the  primary  and  secondary  absorbed  dose  per  incident  parti¬ 
cle  (total  radiation  optional)  for  a  given  spacecraft; 

5)  Determination  of  the  total  radiation  dose,  radiation  dose  rate,  and  LET 
spectrum  for  a  given  space  mission  and  a  given  spacecraft  by  forming  the  product  of 
the  particle  flux  and  the  unit  flux  dose. 

Each  section  in  this  volume  is  self-contained  in  order  to  develop  the  topic  and 
show  its  relationship  to  the  over-all  computer  program  for  calculating  the  radiation 
dose,  dose  rate,  and  LET  spectrum  incurred  in  space  missions. 
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2.0  SPACECRAFT  TRAJECTORIES 


The  space -radiation  flux  is  a  sensitive  function  of  position  and  time;  hence  the 
radiation  flux  encountered  during  a  mission  depends  critical  1/  on  the  accuracy  with 
which  the  spacecraft  trajectory  is  specified.  A  number  of  methods  may  be  employed 
to  calculate  the  spacecraft  trajectory  for  a  given  mission,  depending  on  the  accuracy 
and  detail  of  the  information  required.  Experience  indicates  that  the  required  order 
of  accuracy,  compatible  with  the  variability  of  the  radiation  environment,  is  gen¬ 
erally  ±  tens  of  kilometers  in  altitude,  ±  one  degree  in  latitude  and  longitude,  and 
±  one  minute  in  time;  these  values  are  also  compatible  with  present  uncertainties  of 
flux  values  in  the  radiation  belts. 


g 

An  exact  method,  accurate  to  within  1  part  in  10  but  requiring  much  computer 
time,  can  be  generated  by  numerically  integrating  the  equations  of  motion  (Refer¬ 
ence  1),  considering  the  accelerations  caused  by  the  moon,  the  sun,  and  major 
planets,  the  atmosphere,  and  the  higher  harmonics  of  the  earth's  gravitational  field 
to  an  equal  accuracy.  To  reduce  computing  time,  an  approximate  method  for  tra¬ 
jectory  determination  has  been  adopted  for  radiation -flux  calculations. 


A  second  method,  which  accepts  trajectories  given  point  by  point,  has  been 
included  to  treat  certain  missions,  such  as  powered  flight,  which  cannot  be  deter¬ 
mined  by  the  approximate  method.  In  these  cases,  the  input  trajectories  are  usually 
those  that  have  been  computed  for  other  mission -planning  purposes. 


2.1  APPROXIMATE  TRAJECTORIES 


An  earth  satellite  follows  an  orbit  that  can  be  approximated  by  a  perturbed 
Keplerian  orbit.  The  magnitude  of  these  perturbations  depends  on  the  initial  condi¬ 
tions  of  the  orbit  and  the  time.  The  main  perturbing  Influences  for  near-earth  orbits 
are  the  earth's  oblateness,  the  earth's  atmosphere,  and  the  moon.  These  influences 
cause  periodic  as  well  as  secular  changes  in  the  Keplerian  elements  of  the  orbit.  The 
periodic  changes  can  be  neglected  because  their  amplitudes  are  small;  the  secular 
change  increases  linearly  with  time  and  cannot  be  neglected  for  missions  in  excess  of 
a  few  hours. 


The  effect  of  the  atmosphere  is  to  reduce  the  semimajor  axis  of  the  orbit 
(Reference  2).  In  general  this  effect  Is  bypassed  because  the  radiation  environment 
is  sharply  attenuated  in  the  atmosphere. 

The  moon  exerts  an  average  torque  on  the  orbit  causing  it  to  process  about, 
approximately,  the  pole  of  the  ecliptic.  The  resulting  precession  increases  with  the 
cube  of  the  semimajor  axis  until  it  is  the  dominant,  although  small,  perturbation  out¬ 
side  the  24-hour  orbit.  Since  the  radiation  environment  is  slowly  varying  with 
position  at  these  high  altitudes,  the  effects  caused  by  the  moon  will  be  neglected. 
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Highly  eccentric  orbits,  which  penetrate  close  to  the  earth,  do  experience  a  perigee 
variation  of  the  order  of  a  kilometer  per  day  (Reference  3).  Even  this  variation  is 
not  considered  to  lead  to  sufficient  change  in  the  encountered  flux  to  warrant  Its 
Inclusion. 

The  secular  changes  due  to  the  earth's  oblateness  are  dominant  for  all  orbits 
that  encounter  a  significant  position-dependent  environment.  The  effect  of  the 
higher  harmonics  falls  off  more  rapidly  with  altitude.  Also  the  coefficients  Jj  of  the 
harmonic  terms  with  I  >  2  are  of  the  order  of  the  square  of  J2  (Reference  4).  Since 
the  oblateness  term  yields  a  latitude-longitude  correction  which  is  less  than  1  degree 
per  orbit,  the  higher  order  terms  yield  corrections  of  the  order  of  1  degree  per  300 
orbits.  Accordingly,  all  corrections  except  the  earth's  oblateness  correction  will  be 
neglected  and  this  correction  will  be  treated  only  to  the  first  order  in  J2. 

The  secular  changes  In  the  elements  due  to  the  earth's  oblateness  are  well 
known  (Reference  4).  They  are: 

1)  The  reduction  in  the  period  of  the  orbit; 

2)  The  regression  of  the  Itne  of  nodes; 

3)  The  precession  of  the  line  of  apsides  in  the  orbital  plane. 

To  this  order,  the  geographic  position  and  time  after  Initial  perigee  of  an  earth  satel¬ 
lite  are  given  by: 


r  = 

j 

1  +  e  cos  f 

0  = 

cos  ^  C  sin  i  sin  (9  +  8)  ] 

f  ft  +  tan  ^  [  cos  1  tan  (9  +  B)  ]  , 

cos  (f  +  8)  >  0 

<p  = 

[  ft  +  n  +  tan  ^  E  cos  i  tan  (f  +  6)  3  , 

cos  ( f  +  6)  <  0 

t  «  L  M  =  £  (E  -  e  sin  E) 

Zit  Zir 

where  r,  0  and  9  are  the  geocentric  radius,  colatitude,  and  longitude  of  the 
satellite;^,  l,  «,  and  P  are  the  semilatus  rectum,  the  inclination,  the  eccentricity 
and  the  anomalous  period  of  the  orbit;  ^  M,  and  E  are  the  true,  the  mean,  and  the 
eccentric  anomaly;  ft  is  the  longitude  of  the  ascending  node;  and  8  is  the  argument 
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of  perigee.  The  relationship!  of  these  quantities  (Figure  1)  to  the  input  quantities 
(apogee  altitude  ha/  perigee  altitude  hp,  initial  latitude  Aq,  initial  longitude  $Q) 
are  given  by  (Reference  4): 


ke  -  1.9965  x  10“2(Mm)3//2/sec 


B  -  (30  +  J'ket(R Ji)2  (2  - 1  sin2  i)/a3/2 
J'  -  (3/2)  (Req/Re)2  J2  »  1.63197  x  10“3 
fi  ■  Qo  "  "e^  "  J'keKR ©//)2  cos  i/a3//2 
■  4>0  -  tan  1  (cos  i  tan  Bq) 

Bq  *  sin  ^(cos  9^/sin  i) 

where  the  initial  heading  determines  the  quadrant  of  the  inclination  angle  i.  For 
example,  a  southeasterly  initial  heading  places  the  inclination  angle  in  the  fourth 
quadrant.  The  quantities  R®  and  R«q  are  respectively  the  mean  and  equatorial  radii 
of  the  earth. 

These  expressions  give  the  approximate  position  for  bound  orbits.  The  circular  orbit 
case  follows  from  these  by  setting  e  ■  0.  Since  the  form  is  quite  simple  for  this  case, 
extraneous  computations  of  £,  c,  a,  and  E  are  avoided  by  programming  the  circular 
case  separately. 
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2.2  POINT-BY-POINT  TRAJECTORIES 


Trajectories  that  cannot  be  approximated  by  perturbed  Keplerian  orbits  include 
powered-flight  trajectories,  highly  eccentric  orbits,  and  trajectories  for  missions 
that  have  unique  requirements  of  position  and  time.  Such  trajectories  are  computed 
for  other  mission  requirements,  such  as  power  optimization  and  rendezvous  position, 
and  can  be  used  as  input  for  the  radiation -hazard  trajectory  code.  Thus,  the  trajec¬ 
tory  code  provides  for  the  acceptance  or  generation  of  all  conceived  mission 
trajectories.  To  determine  dose  and  dose  rates,  these  trajectories  are  converted  to 
magnetic  coordinates. 

2.3  CONVERSION  TO  MAGNETIC  COORDINATES 

The  magnetic  coordinates  that  regularize  the  flux  data  best  are  the  magnetic 
flux  density  B  and  the  shell  parameter  L  developed  by  Mcllwain  (Reference  5). 
Alternatively,  the  invariant  radijs  R  and  the  invariant  latitude  X  may  be  useful  for 
simple  geometric  visualization.  For  trapped  radiation  the  validity  of  the  B-L  system 
is  directly  related  to  the  first  and  second  adiabatic  invariants  of  motion  of  the 
trapped  particles.  Therefore,  it  is  assumed  that  the  environment  is  given  in  terms  of 
B  and  L.  In  addition,  it  is  assumed  that  the  magnetic  field  can  be  represented  by  a 
combination  of  n-tuples  of  the  order  6  or  less;  in  particular  the  48  gaussian  coef¬ 
ficients  of  Jensen  and  Cain  (Reference  6).  Mcllwain's  program,  adapted  and  sup¬ 
plied  by  Hess  with  the  Goddard  Space  Flight  Center's  code  designation  C-3,  is  used 
to  calculate  L.  The  method  is  to  calculate  the  integral  Invariant  of  Northrup  and 
Teller  (Reference  7)  and  convert  this  to  the  shell  parameter  L,  with  a  control  error 
ERR.  ERR  is  the  average  fractional  error  in  L  (Reference  8).  The  typical  error  is  0. 1 
ERR,  while  about  i  in  100  cases  will  have  fractional  errors  of  10  ERR.  The  first- 
generation  code  used  ERR  =  0.03.  Since  this  error  can  be  expected  to  be  too  large 
for  some  missions  with  unique  requirements,  the  present  code  has  been  set  to  ERR  * 
0.01.  The  change  added  about  1  minute  of  computer  time  per  trajectory  day  to  the 
existing  code.  Test  runs  indicated  that  roughly  60  percent  of  the  run  time  was  used 
in  calculating  B  and  L.  Therefore,  an  efficient  program  avoids  extraneous  calcula¬ 
tions  of  B  and  L.  A  systematic  study  of  the  B-L  program  is  outside  the  scope  of  this 
work.  However,  as  more  geomagnetic  data  become  available,  such  a  study  should  be 
undertaken.  For  our  purposes  we  assume  that  the  B  and  L  values  can  be  calculated  to 
any  degree  of  accuracy  by  adjusting  ERR  so  that  it  is  consistent  with  the  most  strin¬ 
gent  minion  requirements.  During  the  present  investigation,  where  the  environment 
does  not  warrant  such  accuracy,  methods  were  sought  to  reduce  the  computing  time 
by  avoiding  the  calculation  of  B  and  L  with  the  C-3  code. 

2.4  TIME-SAVING  TECHNIQUES 

The  most  obvious  method  to  decrease  computer  time  is  to  bypass  Mcllwain's 
program  where  the  flux  is  known  to  be  negligible.  For  low -altitude  orbits  the  flux  of 
trapped  radiation  is  zero  for  more  than  half  the  geocentric  spherical  surface.  With  a 
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limited  survey  of  the  low-oltltude  flux  data  (Reference*  9  and  10)  one  deduce*  that 
the  flux  of  trapped  particle*  for  altitude*  under  800  km  i*  confined  to  the  Atlantic 
between  110  degree*  We*t  to  80  degree*  Ea*t.  For  altitude*  under  400  km  the  region 
*hrink*  *omewhat  *o  that  it  is  contained  between  75  degree*  West  and  60  degree* 

East.  Therefore,  if  one  bypasses  the  calculation  outside  this  region,  a  50-percent 
reduction  in  run  time  will  be  realized.  Since  the  shape  and  location  of  this  region 
are  only  approximately  known,  we  make  it  generously  large  and  set  the  flux  of 
trapped  particles  equal  to  zero  for  longitudes  west  of  120  degrees  and  east  of  90 
degrees  East  for  altitudes  below  800  km.  An  option  to  bypass  this  geographic  test  is 
provided  because,  if  one  is  interested  in  the  untrapped  radiation,  magnetic  coordi¬ 
nates  must  be  calculated  at  all  longitudes.  In  regions  suspected  of  having  significant 
radiation,  whether  trapped  or  untrapped,  other  means  of  saving  computer  time  must 
be  employed. 

A  very  promising  method  of  computing  B  and  L  efficiently  is  based  upon  using 
the  C-3  code  only  at  controlled  points  along  the  trajectory.  Let  B.  =  B(tj)  and 
Lj  =  L(tj)  be  the  values  of  the  smooth  functions  B(t)  and  L(t)  at  n  selected  points  tj; 
1*1,  2,  ...  n;  in  an  interval  t]  <  t  <  tn.  Through  the  n  computed  values  of  Lj,  for 
example,  there  is  a  unique  polynomial^  of  order  n-1.  This  polynomial  approximates 
the  shell  parameter  L(t)  in'  the  interval  t]  <  t  <  tp  and  has  been  used  to  calculate 
what  may  be  called  interpolated  values  within  each  interval  At  =  t-+j  -  tj  at  a  minor 
cost  in  computing  time  compared  with  the  time-consuming  integrals  of  the  C-3  code. 
The  magnetic  induction  B(t)  may  be  similarly  treated.  For  a  sample  trajectory. 

Table  1  shows  the  average  percentage  differences  between  the  interpolated  values  at 
seven  points  within  each  Interval  At  and  calculated  values  at  the  same  points  accord¬ 
ing  to  the  C-3  code  with  ERR  *  0.01  for  several  interval  sizes.  The  column  headed 
No.  of  Cases  denotes  the  number  of  intervals  of  size  at  At  considered.  The  averages 
are  then  taken  over  points  which  are,  in  number,  seven  times  the  number  of  cases. 

The  quantities  called 


are  the  maximum  over  the  number  of  cases  of  averages  taken  in  each  case. 
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Table  1 


AVERAGE  AND  MAXIMUM  AVERAGE  PERCENTAGE  DIFFERENCES  OF 
COMPUTED  AND  INTERPOLATED  MAGNETIC  COORDINATES 


At 

(sec) 

No.  of 
Cases 

(percent) 

M 

(percent) 

((■T-),) 

(percent)  m°X 

m,) 

l  l\  max 

(percent) 

60 

7 

.0001 

.0061 

.0006 

.0147 

80 

5 

.0002 

.0036 

.0006 

.0110 

120 

7 

.0017 

.0056 

.0063 

.0203 

160 

10 

.0084 

.0127 

.0236 

.0262 

320 

4 

.0460 

.2050 

.076 

.380 

960 

3 

6.3 

8.6 

12.6 

15.1 

That  ts: 


where  N  denotes  the  number  of  cases  and  B-  denotes  the  computed  B  values  at  the 
k™  point  In  the  j*h  case. 

These  results  indicate  that  the  Interpolation  routine  maintains  the  accuracy  of 
the  C-3  code  for  small  time  intervals.  This  test  run  is  too  limited  to  justify  a  blanket 
recommendation  to  use  this  Interpolation  method  for  time  intervals  less  than  a  speci¬ 
fied  single  value.  However/  for  any  region  of  space,  some  range  of  time  intervals 
exists  where  the  interpolation  method  will  suffice  and  the  C-3  code  may  be  passed. 
Hence,  there  is  promise  of  securing  any  required  density  of  trajectory  points  in  mag¬ 
netic  coordinates  with  little  more  machine  time  than  that  required  by  the  former 
program  (Reference  11).  The  requirements  on  the  density  of  points  is  imposed  by  the 
flux  and  dose  integration,  which  will  be  treated  in  Section  8.0.  A  discussion  of 
various  aspects  of  the  radiation  environment  relating  to  the  dose  rate  will  be  useful 
in  understanding  these  requirements. 


9 


3.0  RADIATION  ENVIRONMENT 


The  ionizing  radiatioh  environment  in  cis-lunar  space  may  be  categorized  as 
geomagnetical  I  y  trapped  or  untrapped  radiation.  In  principle,  the  trapped  radiation 
consists  of  charged  particles  that  remain  in  periodic  oibits  in  the  magnetosphere  for 
several  drift  periods.  In  practice,  however,  this  categorization  is  subjective;  most 
radiation  detectors  do  not  distinguish  between  particles  in  trapped  vrbits  and  those 
that  are  momentarily  trapped.  The  raw  data  obtained  in  the  Van  Allen  belts  is  pro¬ 
cessed  to  obtain  a  flux  which  must  be  considered  trapped  radiation;  that  is,  the  galac¬ 
tic  background  is  subtracted,  the  transient  behavior  is  separated,  and  the  data  are 
smoothed  In  B-L  coordinates.  Thus  the  resulting  flux  is  the  best  present  estimate  of 
the  trapped  component.  All  other  radiation,  such  as  galactic -cosmic  and  solar- 
particle  radiation  will  be  treated  as  untrapped  radiation*  The  level  of  the  ionizing 
neutral  components,  chiefly  photons  and  neutrons,  is  so  low  that  they  are  neglected 
in  this  radiation  hazard  study. 

3.1  TRAPPED  RADIATION 

The  trapped -particle  flux  is  given  experimentally  in  B  and  L,  the  magnetic  co¬ 
ordinates  of  Mcllwain  (Reference  5).  B  is  the  magnetic  induction  at  the  point  in  space. 
L  is  the  magnetic  shell  parameter  that  labels  the  shell  upon  which  the  guiding  center 
of  the  trapped  particles  is  adiabatically  confined  as  it  drifts  around  the  earth. 

Mcllwain  has  shown  that  L  is  approximately  constant  along  the  real  magnetic  field 
line.  In  a  dipole  field,  L  is  rigorously  constant  along  the  field  line  and  has  the  geo¬ 
metric  property  of  being  the  equatorial  distance  from  the  dipole  center  to  the  field 
line.  The  B  surfaces  and  L  surfaces  for  a  dipole  field  are  shown  in  Figure  2.  For  the 
real  field  these  smooth,  symmetrical  surfaces  become  warped  and  asymmetrical  par¬ 
ticularly  at  low  altitudes.  Because  of  the  connection  to  the  adiabatic  invariants  of 
the  trapped -particle  motion,  the  trapped  particle  flux  Is  a  function  of  B  and  L. 

For  analysis,  all  reported  fluxes  are  considered  to  arise  from  particles  adiabatical¬ 
ly  trapped  for  several  drift  periods.  The  characteristics  of  this  flux  that  are  important 
for  dose  calculations  are  the  omnidirectional  flux,  the  spectral  distribution,  and  the 
angular  distribution  at  any  point. 

The  most  common  presentation  of  the  omnidirectional  fluxes,  the  isoflux  contours 
in  B-L  space,  is  called  an  omnidirectional  flux  map.  This  map  is  converted  to  tabular 
form  for  the  computer.  The  simplest  tabulation  to  use  is  that  of  the  fluxes  at  B-L  grid 
points.  Whether  the  conversion  is  accomplished  automatically  or  by  hand,  the  most 
efficient  grid  must  be  determined  from  the  flux  map.  In  general,  this  grid  is  not 
uniform.  A  particularly  successful  method  of  conversion  proceeds  as  follows:  on  each 
of  a  large  number  of  L  shells  the  logarithm  of  the  flux  is  plotted  versus  B.  The  density 
of  B  grid  lines  is  then  chosen  high  enough  to  enable  linear  Interpolation  to  reasonably 
approximate  the  logarithm  of  the  flux  at  Intermediate  B  values  for  all  chosen  L  shells. 
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FIGURE  2.  The  geometry  of  the  B-L  coordinate  system 
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Then  the  logarithm  of  the  flux  at  each  B  grid  line  is  plotted  as  a  function  of  L  and  the 
density  of  L  grid  lines  is  chosen  to  validate  linear  interpolation.  These  plots  are  then 
used  to  extrapolate  in  B  and  L  beyond  the  magnetic  equator  so  as  to  completely  sur¬ 
round  the  magnetic  equator  by  grid  points  with  non-zero  fluxes.  This  ensures  that  the 
linear  interpolation  of  the  logarithm  of  the  flux  will  give  a  reasonably  accurate  value 
of  the  flux  throughout  the  map.  The  flux  maps  are  reported  to  be  accurate  within  a 
factor  of  two.  Since  the  above  conversion  to  tabular  form  can  easily  be  accomplished 
while  maintaining  a  10  to  20  percent  accuracy,  higher -order  interpolation  schemes 
are  not  necessary.  The  proton  map  constructed  by  W.  Hess  from  the  1958  measurements, 
and  the  Aerospace  proton  and  electron  maps  compiled  by  J.  Vette  (Reference  12)  from 
the  1962  measurements  are  used  in  the  program.  Other  tabulations  can  be  substituted 
If  desired.  Since  the  grid  and  extent  of  a  new  flux  table  may  be  different,  such  In¬ 
formation  must  be  provided  as  flux  table  input  data. 

The  angular  distributions  presented  in  this  section  are  based  on  the  maps  of  Hess 
and  Vette;  furthermore,  the  energy  spectra  were  synthesized  to  be  appropriate  for 
these  maps. 

3.1.1  Omnidirectional  Flux  Maps 

The  omnidirectional  flux  is  the  count  rate  of  a  spherical  counter  of  unit  cross  sec¬ 
tional  area  and  of  unit  efficiency.  The  omnidirectional  flux  of  each  particle  type 
measured  in  the  magnetosphere  is  reported  in  B  and  L,  as  determined  by  Mcllwain's 
code.  The  low  altitude  is  averaged  in  longitude;  additionally  the  data  on  high  L- 
shells  are  averaged  in  time. 

The  low-altitude  smoothing  affects  the  fluxes  that  are  at  least  3  orders  of  magni¬ 
tude  down  from  the  maximum  values.  The  longitudinal  variations  occur  In  regions  of 
B  and  L  that  map  into  the  earth  and  dense  atmosphere  in  the  Atlantic.  The  minimum 
altitude  of  a  B-L  trace  is  called  hmfn.  The  B-L  plots  for  a  family  of  hmjn  are  shown 
In  Figure  3.  The  particles  at  B  values  lying  above  the  hmjn  ”  100  curve  will  be 
severely  scattered  over  the  anomaly.  Hence,  they  do  not  remain  trapped  for  a  drift 
period.  The  errors  introduced  by  treating  low-altitude  fluxes  as  trapped  particles 
will  be  insignificant  for  the  normal  trajectory  which  also  passes  through  the  South 
Atlantic  region  where  the  fluxes  are  much  higher. 

The  high  L  shell  smoothing  effects  large-magnitude  fluxes.  The  diurnal  varia- 
tionsr  which  can  amount  to  an  order  of  magnitude,  have  been  averaged  out.  In  most 
instances,  the  periodic  variations  that  can  be  correlated  to  the  sun's  rotation  period 
have  also  been  averaged  out.  These  also  are  an  order  of  magnitude  at  L  “  4  (Refer¬ 
ence  13).  Although  such  data  give  only  the  average  flux  and  dose  rate,  these 
averages  are  useful  estimates  in  any  case  and  become  correct  for  long-term  orbital 
missions.  A  discussion  of  the  time  variations  in  flux  and  dose  rate  which  chiefly  af¬ 
fect  short  missions  is  given  In  Section  3. 1.3. 
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L  (earth  radii) 

FIGURE  X  A  B-L  plot  of  the  minimum  altitude  h  . 


The  fl  ux  map  compiled  by  Hess  from  data  obtained  by  Explorer  IV  (Reference  5) 
and  Pioneer  III  end  IV  is  shown  in  Figure  4.  The  contours  are  labeled  by  the  integral 
omnidirectional  flux  of  protons  of  energies  greater  than  30  Mev.  The  composite  spec¬ 
trum  normalized  to  1  above  30  Mev  which  was  proposed  by  Hess  is  given  by: 


j(E)-  0.009019 

-  0. 01443 lexp(-E/85) 

-  0. 0044095exp  (-E/172) 
■=  0. 0022047exp  (-E/245) 


20  <  E  <  40  Mev 

40  <  E  <  200  Mev 
200  <  E  <  400  Mev 
400  <  E  <  1000  Mev 


This  spectrum  is  portrayed  graphically  in  Figure  5.  The  high-energy  proton  map  com¬ 
piled  by  Vette  (Reference  12)  from  data  obtained  by  Mcllwain  on  Explorer  XV  aug¬ 
mented  by  data  of  Freden  (Reference  10)  is  shown  in  Figure  6.  The  contours  are 
labeled  by  the  integral  flux  from  40  to  80  Mev.  The  spectrum  proposed  by  Vette 
normalized  to  one  between  40-80  Mev  is  given  as: 


j(E)»  0.0314  20<E<  40  Mev 

-  0.  Q505expf  E/84)  40  <  E  <  1 50  Mev 

-  0. 0276exp  (-E/1 25)  1 50  <  E 

This  spectrum  is  compared  with  Hess'  in  Figure  5.  The  electron  flux  map  representing 
the  environment  on  about  1  November  1962  was  compiled  by  Vette  and  is  portrayed 
in  Figure  7.  The  contours  are  labeled  by  the  integral  flux  of  electrons  of  energies 
greater  than  0.  5  Mev.  The  spectrum,  also  proposed  by  Vette,  normalized  to  one 
above  0.  5  Mev  Is  given  as: 

j(E)  «  1. 024exp  (-0. 575 E  -  0. 055E2)  L  <  1. 5 

■  4. 082  exp  (-E/0. 58)  L  >  1.5 

This  spectrum  is  portrayed  in  Figure  8.  The  low-energy-proton  map  compiled  by  Vette, 
based  on  Freden's  data  in  1962,  is  shown  in  Figure  9.  The  contours  are  labeled  by  the 
integral  flux  of  protons  of  energies  between  5  and  20  Mev.  The  spectrum  that  Is 
suggested  by  Vette,  normalized  to  one  between  5  and  20  Mev  is  given  by: 

j(E)  -  974. 79 E"4, 5  4  <  E  <  15  Mev 

-  974. 79  (15f 4* 5  -  4. 972  x  10"3  15  <  E  <  20  Mev 
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L  (earth  radii) 

FIGURE  4.  The  Hess  P-1  proton  flux  map 


Protons  (cm  sec  Mev 


L  (earth  radii) 

FIGURE  6.  The  Vette  high-energy  proton  map 


0.40 
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E  (Mow) 


FIGURE  8.  The  electron  spectrum  proposed  by  Vette 
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L  (earth  radii) 

FIGURE  9.  The  Vette  low-energy  proton  mop 


Use  of  these  spectra  would  yield  a  one-to-one  correspondence  between  the  om¬ 
nidirectional  flux  and  the  dose  rate  Inside  a  uniformly  shielded  chamber.  Dosimetry 
measurements  (Reference  14)  show  that  there  Is  no  frne-to-one  correspondence  and 
that  spectral  and  angular  variations  must  be  considered. 

3.1.2  Spectral  Dependence 

Direct  spectral  measurements  of  the  protons  (Reference  15)  and  of  the  electrons 
(Reference  16)  have  established  that  the  spectrum  Is  strongly  L-dependent  and  at 
least  slightly  B-dependent.  Additional  analysis  of  Explorer  IV  data  (Reference  17) 
and  of  Midas  data  (Reference  18)  supports  this  conclusion.  To  determine  the  dose 
rates  at  any  point  It  is  necessary  to  know  the  appropriate  energy  spectrum  for  each 
particle  type  as  well  as  its  respective  flux.  Hence,  a  radiation  environment 
program  must  specify  the  spectra  as  well  as  the  fluxes^of  all  particles  at  a  particular 
B-L  point.  In  addition,  this  specification  must  be  In  a  form  In  which  the  flux  can  be 
Integrated  from  point  to  point,  retaining  the  spectral  dependence  so  that  the  total 
differential  energy  flux  for  any  mission  can  be  determined. 

It  is  possible,  of  course,  to  specify  the  differential  energy  flux  at  any  point  on 
an  energy  grid  and  then  simply  integrate.  In  mission  time,  these  energy-grid  fluxes. 

To  calculate  the  dose  rate  and  mission  dose  from  these  data  would  require  that  the 
many  spectra  generated  be  processed  by  the  shielding  code.  This  extremely  costly 
process  requires,  at  the  very  minimum,  20  separate  runs,  one  for  each  grid  point.  It 
appears  more  efficient  to  fit  the  measured  spectrum  at  each  point  to  a  linear  combina¬ 
tion  of  a  finite  number  (<  10)  of  basic  spectra  and  determine  the  flux  and  total  mis¬ 
sion  flux  of  each  spectral  type.  Then  the  dose  rate  and  mission  dose  are  easily 
obtained  through  multiplication  by  the  respective  dose  conversion  numbers  obtained 
with  a  finite  number  (<  10)  of  runs  of  the  shielding  code  and  through  summing  these 
contributions.  In  addition,  if  the  basic  spectra  are  chosen  to  be  characteristic  of 
the  sources  and  of  the  predominant  modulations  by  the  loss  mechanisms,  this  decom¬ 
position  of  the  spectral  dependence  might  aid  attempts  at  understanding  these 
mechanisms.  Unfortunately,  the  source  and  loss  mechanisms  have  not  been  well 
enough  determined  to  realize  this  additional  value. 

The  electron  spectra,  as  measured  by  H.  West  (Reference  16)  in  the  late  fall  of 
1962,  provide  a  fair  coverage  of  the  electron  belt.  The  measurements  were  performed 
with  a  five-channel,  permanent -magnet,  180-degree -focusing  beta-ray  spectrometer. 
The  five  energy  channels  were  0.325,  0.955,  1.63,  2.40  and  3.25  Mev.  Sixty  rep¬ 
resentative  spectra,  measured  in  the  range  B  <  0. 24  gauss  for  L  <  1.7  and  B  ~  0.010 
gauss  for  L  >  1.7,  were  kindly  supplied  by  H.  West.  These  measurements  indicated 
that  the  electron  spectrum  softens  from  a  fission-like  spectrum  as  L  increases  out  to 
about  1.6.  Here  an  Inflection  in  the  spectrum  between  1  and  2  Mev  is  apparent  and 
is  portrayed  as  a  depletion  In  the  flux  near  1  Mev  as  L  increases  to  2.7.  As  L  in¬ 
creases  from  2.7,  the  flux  at  1  Mev  builds  back  up  while  the  spectrum  above  2.4  Mev 
becomes  steeper  until  at  L  ■  4  the  spectrum  resembles  that  at  L  =  1.5.  Then 
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as  L  increases  above  4,  the  spectrum  continually  softens.  At  low  L  values  and  large 
B  values  the  differential  flux  was  depleted  at  low  energies  in  a  manner  consistent 
with  the  calculations  of  Walt  (Reference  19)  on  the  atmospheric  scattering  loss 
mechanism. 


On  the  basis  of  this  qualitative  behavior/  the  spectre  measured  at  L  >  1.7  on 
October  27,  1962/  were  processed  utilizing  five  basic  soectra.  These  measured 
spectra  contain  an  undetermined  contribution  from  a  Russian  high -altitude  test  con¬ 
ducted  earlier  in  October.  These  five  spectra  included  as  source-characterizing 
spectra  the  Carter  (Reference  20)  fission  spectrum/  a  steep  low-energy  spectrum  ex¬ 
pected  from  a  neutron  albedo  source,  end  an  exponential  spectrum  with  a  mean  energy 
of  0.5  Mev  for  the  outer  belt.  The  loss-modulated  spectra  were  represented  by  a 
Gaussian  idealization  of  the  spectrum  measured  by  West  at  L  *  1. 15,  B  *  0.258,  and 
a  spectrum  resembling  the  one  measured  at  L  M  2. 7,  B  ■  0.08.  Linear  combinations 
of  these  five  were  least-squares  fitted  to  the  smoothed  spectra  of  West  at  every  0.5 
Mev  from  0.5  to  4.0  Mev.  Although  the  fit  in  this  range  was  promising  (relative 
standard  deviations  a  on  the  eight  points  from  0.01  to  0.2),  the  spectrum  was  unreal¬ 
istic  outside  this  range,  in  a  number  of  cases.  This  can  be  expected;  even  a  negative 
differential  flux  can  occur.  For  example,  suppose  the  best  fit  in  the  1-  to  3-Mev 
range  requires  a  negative  coefficient  for  the  fission  spectrum;  then,  since  the  fission 
spectrum  is  the  hardest  spectrum  of  the  set,  the  resultant  fitted  spectrum  must  become 
negative  at  some  higher  energy. 

Another  disturbing  aspect  of  the  use  of  this  set  occurred  when  a  fit  to  the  spectra 
measured  at  L  values  less  than  1.7  was  attempted.  Although  the  fit  was  reasonable 
(o  <  0.25),  none  of  the  spectra  were  as  well  fit  as  in  the  slot  region,  2.0  <  L  <3, 
where  a  *  0.01.  Evidently,  this  basic  set  was  securing  the  best  fit  where  the  flux 
and  the  expected  dose  rate  were  relatively  small. 


The  complexity  of  the  spectral  variations  at  low  L  values,  which  are  evident  in 
the  60  representative  experimental  spectre  plus  the  promise  of  better  definition  at  low 
altitudes  as  more  of  West's  data  are  reduced,  called  for  a  more  flexible  basic  set. 

The  following  set  of  eight  spectre  ^normalized  to  unity  above  0.5  Mev,  were  chosen: 


f-j(E) 


39.745 


r  exp  (-7.35E) 

6.093  x  105exp(-24E) 


E  <  0. 8  Mev 
E  >  0.8  Mev 


f2(E)  -  5.2603exp[-5.555(E  -0.3)23 
f3(E)  -  1.4633exp  t-5.555(E  -0.9)21 
f4(E)  -  1.3303exp  t-5.555(E  -  1.5)2  3 
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f5(E)  -  1.3298exp[-5.555(E  -2. 1)2 3 
f6(E)  -  1.3298exp[-5.555(E  -2.7)2] 

fexp[-5.555(E-3«3)23  E  <  3.3  Mev 
3, 827. 5exp  (-2. 5E)  E  >  3. 3  Mev 


f  exp  [-5. 555(E  -  3. 3)2  3  E  <  3. 3  Mev 

ffl(E)  -  0.77189 

0  1 2. 1 39exp  (bOLSZS  E«-  0. 0055EZ)  E  >  3. 3  Mev 

These  eight  spectra  are  displayed  in  Figure  10. 

To  reduce  the  probability  of  obtaining  unrealistic  differential  flux  at  high  and 
low  energies,  W est*s  spectra  were  extrapolated  down  to  zero  energy  and  out  to  7  Mev. 
A  15-point  least-squares  fit  (at  each  0.5  Mev)  was  then  obtained  for  the  amplitudes 
of  these  eight  basic  spectra.  Thus  each  of  the  60  spectra  Gj(E)  Is  fitted  to  an  an¬ 
alytic  function  of  the  forms 


is  minimized.  Nearly  all  values  of*j  were  in  the  range  0.01  <  <  0.04. 


Each  spectrum  G.(E)  was  based  on  the  measured  spectrum  at  a  point  (B  ,  L  ). 

The  resulting  coefficients  Cj(()  were  normalized  to  a  unit  sum.  That  is,  coefficients 
Cj(B  ,  L  )  were  defined  such  thats 


Hence: 


ci<VLi> 


These  coefficients,  grouped  into  the  B  ranges  0  to  0. 10,  0. 125  to  0. 17,  0. 17  to 
0.20,  were  plotted  against  L  Smooth  curves  through  the  data  were  then  interpreted 
as  the  L  dependence  of  the  coefficients  on  the  B  grid  lines  Bj  ■  0.08,  0.  14,  and 
0. 18  gauss.  The  density  of  L— grid  lines  was  chosen  so  as  to  validate  linear  inter¬ 
polation.  The  L  dependence  at  the  equator  is  shown  graphically  in  Figures  11  and 
12.  The  B  dependence  is  leu  pronounced  as  seen  from  the  example  in  Figure  13. 

The  resultant  C  tables  of  the  coefficients  C;  (B,  L)  were  then  checked  for  unit 
normalization  and  submitted  as  data  to  the  trajectory  and  environment  code. 

This  detailed  approach  employed  for  electrons  is  not  feasible  for  protons  because 
the  spectral  data  are  so  limited.  Where  available,  measured  spectra  were  chosen  as 
basic  spectral  forms.  On  the  higher  L  shells,  the  direct  spectral  measurements  were 
supplemented  with  some  geiger  counter  data,  in  general,  the  spectrum  of  the  trapped 
protons  softens  with  increasing  L.  Mcllwain  (Reference  17)  deduces  that  the  data 
above  L  ■  1.4  are  consistent  with  a  spectrum  for  protons  above  30  Mev  of  the  form: 


j(E)  -  J-  expt-<E  -  30)/Eo(L)] 
with  Eq  ■  a  A 

A  least  squares  analysis  of  Explorer  IV  data  from  31  to  43  Mev  gave 
a  -  306  ±  28 
b  -  5. 2  ±  0. 2 

Imhof  and  Smith  (Reference  18)  find  that  their  data  are  consistent  with  the  same  form 
If 


a  -  460  ±  110 
b  -  4.8  ±  0.4 

These  are  Geiger  counter  and  scintillation  counter  measurements  and  as  such  are  not 
comprehensive  measurements  of  the  energy  spectrum.  The  only  comprehensive  mea¬ 
surements  of  the  positional  dependence  of  the  proton  energy  spectrum  are  reported 
by  Naugle  and  Kniffen  (Reference  15). 
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As  a  reasonable  description  of  the  spectral  dependence  of  protons  greater  than 
30  Mev,  the  following  spectra  are  proposed. 

1)  For  the  Hess  P-1  Map  (see  Figure  4)  the  following  spectra  on  an  L  grid  are 
to  be  normalized  to  one  above  30  Mev. 

a)  For  L  <1.2  the  composite  spectrum  of  Hess  (see  Figure  5). 

b)  At  L  *  1.3  the  measured  spectrum  of  Freden  and  White  (Reference  21) 
labeled  f2» 

c)  At  L  “  1.47,  1.54,  1.64,  and  1.72  the  measured  spectra  of  Nangle 
and  Kniffen  (Reference  15)  labeled  f3,  f^,  f^  and  f^,  respectively. 

d)  And 

f7H  -  0. 22408exp  (-E/20)  at  L  -  1.9 

f8H  -  00. 645exp(-E/5)  for  L  >  2.5 

The  measured  spectra,  f2  through  f^,  are  given  in  Tabel  2  along  with  the  fac¬ 
tors  fj^/fj  which  normalize  each  to  one  above  30  Mev.  These  spectra  are  also  dis¬ 
played  in  Figure  14.  At  L  values  between  two  grid  values  a  linear  combination  of  the 
grid  spectra  is  used.  For  example,  at  L  “  1.6  the  spectrum  would  be  given  byi 

fH(E)  -  0.4  f4H  +0.6  f5H 

2)  For  the  high-energy  proton  map  (see  Figure  6)  the  following  spectra  on  an 
L  grid  are  to  be  normal  ized  to  one  between  40  and  80  Mev. 

a)  For  L  <  1.2  the  composite  spectrum  of  Vette  (see  Figure  5) 

b)  At  L  ■  1.3,  1.47,  1.54,  1.64,  and  1.72  the  measured  spectra  fjj, 

f-4/  f^  and  fg ,  repsectively. 

c)  And 

f7V  -  0.42727exp(-E/20)  at  L  -  2.0 

fg^  *  597.8exp  (-E/5)  for  L  >  2.5 

V 

The  factors  fj  /fj  which  normalize  the  measured  spectra  to  one  between  40  and  80 
Mev  are  included  in  Table  2.  The  coefficients  C;  are  linearly  interpolated  as  before. 
Thus  at  L  ■  1.6,  C4  *  0.4,  C5  ■  0.6  and  the  remaining  Cj  ■  0. 
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Protons  (cm  sec  1  Mev 
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FIGUtE  14,  The  basic  proton  spectra  used  In  the  program 
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For  the  low-energy-proton  map  (Figure  9)  no  spectral  variation  with  position 
is  given.  Hence  the  recommended  spectrum  for  this  map  is  that  one  given  in 
Section  3. 1.  1.  Therefore  the  C  tables  for  the  proton  maps  will  consist  of  ones  and 
zeros. 


As  more  spectral  data  is  reduced  and  obtained,  the  C  tables  for  the  electron 
and  proton  maps  can  be  improved.  It  is  believed,  however,  that  the  predominant 
spatial  variations  of  the  spectra  have  already  been  incorporated. 

3. 1.3  Time  Dependence 

The  time  dependence  of  the  trapped-radiation  flux  is  not  well  determined. 
Hence  it  is  difficult  and  perhaps  misleading  to  develop  a  specific  subroutine  to  treat 
this  dependence.  Another  barrier  to  such  a  program  involves  a  lack  of  agreement 
as  to  what  types  of  time  dependence  are  most  important  to  a  radiation-hazard  study. 

A  brief  discussion  is  given  of  the  features  of  the  time  variation  of  the  belts  and  a 
model  is  proposed  which  treats  the  predominant  features. 

Theoretically,  the  temporal  variations  of  the  trapped-particle  flux  can  be 
qualitatively  decomposed  into  transient  behavior  and  non-transient  behavior.  The 
transient  behavior  is  expected  following  impulsive  injection  such  as  that  experienced 
during  a  high-altitude  nuclear  detonation  or  following  an  eruption  on  the  sun.  The 
nontransient  behavior  is  expected  from  the  quasi-periodic  variations  of  the  sources  or 
loss  mechanisms  for  the  belts.  In  addition,  an  apparent  nontransient  phenomena  is 
expected  from  the  slowly  Varying  B-L  coordinates  at  a  particular  geographic  point. 

No  quantitative  theory  of  these  variations  has  been  developed. 

Experimentally,  a  separation  of  these  variations  is  hampered  by  the  lack  of  a 
quantitative  theory.  However,  for  L  <  2.5  it  appears  that  the  protons  above  30  Mev 
energy  do  not  experience  any  transient  behavior.  These  protons,  in  comparison  with 
less  energetic  protons,  exhibit  only  a  relatively  slight  solar-cycle  variation.  The 
electron  belts  do  experience  a  complex  time  history.  Van  Allen  (Reference  22), 
Brown  (Reference  13),  and  Pizzella  (Reference  24)  have  observed  the  transient 
behavior  of  the  electron  flux  following  July  9,  1962  —  the  date  of  the  high-altitude 
test  in  the  Pacific— as  well  as  during  geomagnetic  storm  activity  in  the  outer  belt. 

It  was  formerly  supposed  that  the  predominant  variation  would  be  in  the  integral  flux. 
West's  spectral  data  indicate  that  for  L  values  above  1.7  the  spectral  variations  with 
time  are  quite  severe  at  about  1  Mev.  Also  the  low-altitude  data  is  consistent  with 
the  theory  of  Walt  (Reference  19)  which  predicts  time-dependent  spectra.  These 
considerations  invalidate  the  use  of  a  spectrum  that  is  constant.  Pizzella  (Reference 
24)  observed  diurnal  variations  in  the  outer  belt.  Brown  (Reference  23)  has  noted 
quasi-periodic  variations  at  L  *  4  which  have  roughly  a  27-day  period.  Cansidenble 
work  remains  to  determine  the  specific  causes,  although  they  are  evidently  associated 
with  the  sun's  influence. 
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The  importance  of  temporal  variations  in  the  environment  to  a  space-radiation- 
hazard  study  is  intimately  dependent  on  the  duration  of  the  mission  and  specific 
requirements  of  the  mission  upon  the  time-dependent  dose.  Thus,  orbital  missions, 
longer  than  a  day  without  specific  requirements  on  the  maximum  dose  rate  do  not  re¬ 
quire  a  detailed  treatment  of  the  diurnal  variation.  To  evaluate  the  hazard  far 
short-term  missions  and  those  with  specific  requirements  a  program  to  systematically 
evaluate  the  data  on  time  dependence  is  required.  For  example,  the  diurnal  varia¬ 
tion  should  be  decomposed  into  variations  that  are  actually  associated  with  source- 
and-loss  variations  and  those  that  are  associated  with  local-time  variations  in  B  and 
L.  At  this  point  a  code  to  calculate  B  and  L  as  a  function  of  local  time  would  be 
useful. 

At  present,  the  only  reasonable  requirement  on  the  electron  code  is  compati¬ 
bility  with  a  time-dependent  subroutine.  Construction  of  the  explicit  form  of  the 
subroutine  will  have  to  wait  for  better  understanding  of  the  time  dependence. 

Data  on  the  protons  has  been  gathered  for  a  longer  period  of  time.  A  reasonable 
method  of  evaluating  the  solar-cycle  variation  of  the  hazard  from  these  protons  is  as 
follows.  The  flux  map  of  Hess  is  interpreted  as  the  environment  at  solar  maximum 
while  that  of  Vette  is  interpreted  as  the  environment  at  solar  minimum.  The  dose 
derived  for  any  mission  from  Hess'  map  can  be  labeled  D  and  that  derived  from 
Vette's  map  can  be  labeled  D  .  .  These  will  represent  ffie* range  of  the  proton  dose 
for  the  given  mission.  Hopefully,  the  dose  during  the  cycle  can  be  approximated  by 

n/.\  _  max  mm  .  max  min  f  ~ 

D(t)» - j -  cosl  2ir 

For  the  transient  belt  of  artificially  produced  electrons  the  environment  code  is  com¬ 
patible  with  a  subroutine  which  would  modulate  In  time  the  flux  map  or  the  C  tables. 
That  is, the  subroutine  would  calculate  a  time  dependence  of  the  forms 

F(B,L,t)-GBA(t)F(B,l) 

and/or 

Cj(B,L,t).H^  I,  L(t)C.(B,L) 

The  time  modulating  operators  G(t)  and  H^(t)  are  still  to  be  determined  from  the  data. 
Their  determination  should  be  one  of  the  objectives  of  future  work  planning. 

3. 1.4  Angular  Distributions 

The  angular  distribution  of  geomagnetical ly  trapped  particles  is  generally  not 
isotropic,  although  in  some  cases  isotropy  is  a  reasonable  approximation  as  the  follow¬ 
ing  discussion  will  show.  The  following  derivation  of  the  relationships  for  determining 


t  -  t  i 
max 

1 1  years^ 
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the  angular  dependence  of  trapped  radiation  was  undertaken  to  determine  the  effect 
of  angular  dependence  on  the  radiation  dose  within  a  spacecraft.  The  effect  for  a 
particular  case  is  given  in  Section  7.0.  Another  purpose  in  developing  this  code  is 
to  facilitate  the  interpretation  of  experimental  data  and  the  design  of  space-borne 
radiation  detection  systems. 

In  the  steady  state  the  angular  distributioi.  at  any  point  has  the  general  property, 
which  greatly  simplifies  the  description,  that  it  may  be  specified  by  the  distribution 
in  the  pitch  ang[|  a,  the  angle  between  the^velocity  vector  and  the  local  magnetic- 
induction  field  B.  If  the  direction  of  the  B  field  is  considered  as  a  polar  axis,  then 
the  distribution  in  the  azimuthal  angle  is  uniform  (two-dimensional ly  isotropic  ?n  the 
plane  perpendicular  to  the  field  line).  At  low  altitudes  for  particles  of  large  radius 
of  gyration,  the  variation  in  the  atmospheric  density  is  appreciable  over  one  radius  of 
gyration  and  thus  produces  the  "East-West  asymmetry"  (Reference  25).  In  this  treat¬ 
ment  of  angular  dependence  these  atmospheric  effects  will  not  be  considered  nor  will 
a  description  of  the  "edges"  of  the  trapped  belts  be  obtained.  Rather,  effort  will  be 
concentrated  on  the  regions  where  the  adiabatic  invariants  provide  a  working  des¬ 
cription  of  the  particle  motion.  Most  of  the  flux  lies  in  such  regions.  There  the  flux 
is  uniformly  distributed  in  azimuthal  angle  because  the  particles  spiral  about  the  field 
lines  and  the  phase  or  azimuthal  angle  along  this  trajectory  Is  random. 

The  first  part  of  the  discussion  will  show  that  the  pitch-angle  distribution  at  any 
point  along  a  given  field  line,  labeled  by  L,  can  be  determined  from  knowledge  of 
the  flux  as  a  function  of  position  (References  26,  27,  and  28).  The  flux  maps  are 
designed  to  provide  the  latter  quantity.  Thus  the  problems  involved  in  determining 
the  pitch-angle  distributions  were  chiefly  numerical.  They  have  been  solved 
successfully.  Figure  15  Illustrates  the  geometry  and  the  parameters. 

However,  to  evaluate  the  relative  merits  of  various  shielding  configurations 
with  accurate  account  of  the  angular  distributions,  the  distributions  must  be  trans¬ 
formed  at  each  point  into  coordinates  fixed  in  the  space  vehicle.  The  additional  vari¬ 
ables  required  to  keep  track  of  the  orientation  of  the  vehicle  at  all  times  and  to  trans¬ 
form  the  angular  distribution  to  this  coordinate  system  as  a  function  of  time  demands, 
in  general,  an  impractical ly  large  amount  of  time  and  storage  space  in  the  computer. 

Hence,  the  approach  has  been  that  the  routine  will  generate  the  pitch-cngle 
distributions  at  any  or  all  points.  The  comparison  of  shielding  configurations  will 
then  be  made  point  by  point,  by  hand.  Two  somewhat  realistic  mission  plans  can  be 
treated  exactly  by  just  summing  such  data;  one  where  the  vehicle  has  a  fixed  orien¬ 
tation  with  respect  to  the  TT  field,  and  the  otherwhere  the  vehicle  spins  rapidly  about 
an  axis  that  maintains  a  constant  angle  with  the  B  field.  Such  missions  will  usually 
tnclude  the  1»est"  and  "wont"  cases  for  a  given  shield  became  the  angular  distribu¬ 
tions  quite  consistently  are  largest  perpendicular  to  the  field  line  and  decrease  mono- 
tonlcally  towards  the  field  line.  For  example,  some  typical  distributions  at  the 
equator  are  shown  in  Figures  16  through  20.  The  equatorial  distributions  are  the 
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Electrons  /(cm  sec  cos  a)  x  Electrons/(cm  sec  cos  a 


Electrons  Electrons 

FIGURE  20.  Equatorial  pitch  angle  distributions  at  L  «  3. 5  and  l  =  5.  0  for  the  electron  and  low-energy 
proton  maps 


broadest  and  have  the  greatest  chance  to  peak  away  from  perpendicular.  But  none  of 
the  natural  ones  appear  to  have  statistically  significant  peaks  at  angles  other  than  90 
degrees  with  the  field  line. 

The  equations  to  be  derived  here  relate  an  experimentally  measured  quantity, 
the  omnidirectional  flux  F{B,  L)  as  a  function  of  the  B-L  coordinates  to  an  experimen¬ 
tally  measurable  directional  flux  Fa(B,  L).  The  fluxes  are  equal  to  the  products  of 
the  density  and  the  average  speed  (magnitude  of  velocity)  of  the  particles  accepted 
by  each  flux  measurement,  respectively.  If  the  energy  spectrum  is  constant  along  a 
field  line,  then  the  average  speed  is  the  same  for  the  directional  fluxes  and  the 
omnidirectional  flux.  This  is  approximately  true  in  the  belts  described  by  the  flux 
maps.  Thus  there  is  just  a  constant  factor  between  density  and  flux  so  that  in  the 
following  deviations,  density  may  be  replaced  everywhere  by  flux. 

Under  steady-state  conditions,  let  N  be_£he  total  number  of  particles  contained 
in  a  small  flux  tube  labeled  by  L  Since  V*  B  3  0,  the  cross^ectional  area  A  of 
this  tube  will  vary  with  magnetic  induction  as 

BA  -  B0A0  (1) 

where  the  subscript  zero  refers  to  quantities  evaluated  at  the  equator  or  mid-plane. 

A  is  chosen  small  enough  that  the  particle  density  is  uniform  over  A. 

Let  aQ  denote  the  pitch  angle  at  the  equator  and  a  the  pitch  angle  at  some 
other  point  which  is  located  by  its  magnetic  induction.  See  Figure  15  fora  geomet¬ 
rical  Interpretation  of  these  quantities.  Then  conservation  of  magnetic  moment,  an 
adiabatic  approximation,  implies  that  a  particle  follows  a  trajectory  such  that 


sin^ 


sin  a  =  constant 


(2) 


Let  f(p)dp  ■  the  fraction  of  the  N  particles  which  have  their  "magnetic  moment" 
between  p  and  p  +  dp.  The  quantity  p,  called  here  the  magnetic  moment,  is 
really  a  dimensionless  quantity  proportional  to  the  conventional  magnetic  moment. 
The  latter  is  given  by  p  pv/(2  B^,  or  nonrelativistically  by  p  m  v  /(2  B0)  (p  - 
momentum,  m  ■  maw,  v  *  velocity).  The  dimensional  factors  are  constants. 

If  a  is  related  to  p  by  Equation  2,  and  If  |  sin  a  |  <  1  so  that  particles  of 
"magnetic  moment"  p  penetrate  at  least  as  far  as  B  before  mirroring,  then  this 
fraction  f(p)dp  is  equal  to  the  fraction  of  the  N  particles  which  have  the  cosine  of 
their  pitch  angles  between  cosa  and  cosa  +  d(cosa)  when  they  crow  a  section  normal 
to  the  tube  at  B.  Denote  this  latter  fraction  by  fg  (cosa)  d(cosa).  The  previous 
statement  corresponds  to  the  equation 
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(3) 


fg  (cosa)J  d(cosa)|  ■  f(p)|dp| 

Equation  3  really  serves  to  define  fg(cosa). 

The  function  of  primary  interest,  however,  is  a  density  denoted  by  ng(cosa) 
such  that 

ng(cosa)  d(cosa)  dV  ■  the  number  of  Darticles  with  the  cosine  of 

their  pitch  angle  between  cos  a  and  cos  a  + 
d(cosa)  and  in  the  small  volume  dV  about  a 
point  B  at  any  instant  of  time. 

Let 

dV  =  dAds 

where 

dA  ■  small  element  of  A 

ds  =  small  element  of  length  along  the  field  line  at  B 

Then 

-jq-  ng(cosa)  d(cosa)  dV  «  fraction  of  all  particles  in  this  flux  tube  with 

cosine  of  pitch  angle  between  cos  a  and 
cosa  +  d(cosa)  and  which  also  lie  in  the 
volume  dA  ds  at  any  one  timet. 

The  fraction  of  the  particles  with  pitch  angle  a  which  actually  is  In  the  length  ds  of 
the  flux  tube  at  any  one  time  is  just  the  fraction  of  their  trajectory  which  is  in  ds; 
namely, 

dj  m  1  ds 

B)  /(a,B)  cosa  ™ 

where 

Jt  (  a,  B)  -  length  of  the  trajectory  of  particle  at  B  with  pitch  angle  a 
dlm  length  of  the  part  of  the  trajectory  in  ds 
since  the  angle  between  the  trajectory  and  the  field  line  is  simply  a. 
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The  trajectory  length  Is  given  by 


./(a,!*)  =  /— — 

J  cos  a 

but  this  quantity  will  disappear  in  the  final  result.  Actually,  J ox\\y  depends  on  one 
parameter,  which  can  be  taken  to  be  p  ,  which  is  in  turn  a  function  of  a  and  B  by 
Equation  2.  Therefore 

/(a,  B)  =  /(p(a,B)) 

Finally,  since  the  density  is  uniform  across  the  area  A,  the  fraction  of  particles 
in  the  volume  Ads  of  the  tube  and  also  in  dV  =  dAds  is  just  dA/A.  Thus, 

ng(cosa)  d(cosa)  dV  *  (total  number  of  particles)  x 

(fraction  of  all  particles  with  cosine  of  pitch  angle 
between  cos  a  and  cos  a  +  d(cosa))  x 
(fraction  of  these  in  ds  at  any  one  time)  x 
(fraction  of  preceding  fraction  in  dA) 

-  N fg(co.a)  d(co»a)i([j(o>di))  cotn  ^ 


Or 


N  fg(cos  a) 

ng(cosa)  ^|(|1(a/  B))cosa 


(5) 


The  total  number  density  at  B  (an  experimentally  measurable  quantity  —  number 
per  unit  volume)  is  then 


n(B) 


L 


cos  a  L 

ng(cosa)  d(cosa) 


cosol 


fg(cosa)  d(cosa) 
Jf!p(a,6))cosa 


(6) 


where  is  the  minimum  pitch  angle  for  trapped  particles  at  a  given  B  and  L.  The 
integrand  can  be  related  to  the  function  f(p)  by  changing  the  variable  of  integration 
to  p  where  a  and  p  are  again  related  by  Equation  2.  Thus 


r  b 

cos  a  ■  Y  1  -  g—  p 


(7) 
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Equation  3  may  be  used  to  eliminate  the  necessity  for  computing 
That  is. 


d(cos  a) 
dp 


fg(cosa)  |d(cosa)|  *  fg(cosa(p))  |  dp! 

■  M  |  dp  | 

Since  the  limits  of  the  integral  Equation  6  are 


0  <  cos  a  <  cosa^ 


t 


and 


T 


.  2 
sin  a 


-jp(1  -  cos2  a) 


the  limits  on  p  are 


>  M  >  ~g~ 


.  2 

sin  aL 


ML 


(8) 


The  minimum  magnetic  moment  for  trapping  is  independent  of  position  and  hence  pl 
does  not  depend  on  B.  Thus,  substitution  of  Equations  1,  3,  7,  and  8  into  Equation  6 
gives 


where  x(m)  *  f(p )/X  (p)  Is  the  unknown  function  we  wish  to  determine,  because  the 
flux  per  steradian  at  B  with  pitch  angle  a  is 
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Equation  9  is  an  integral  equation  for  x(p)  'n  terms  of  the  measured  density  n(B). 
The  equation  is  of  a  type  called  Abel's  equation  (Reference  29).  To  reduce  (9)  to 
standard  form  let  Bq/B  =  x  and/lT"  nCB^/x)  =  v(x).  Thus 

x 

x(f')  dp* 

\/  XV 

-1/2 

To  solve  for  x  multiply  both  sides  of  the  equation  by  (p-x)  and  integrate  over 
x  from  p|_  t°  F«  Thus 


dx  dp1 


Now  differentiate  both  sides  with  respect  to  p.  Let  dv(x)/dx  =  v'(x). 


Then 


v(x)  dx 


N  J_ 
9p 


|  h  rv 

-  2^/p -x)x  n(Bo/x)  +  2  /  v'(x)  dx  | 

*FL  'ML 
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Since  n(B<y/pL) 

x'(m) 


y  f  >/»-*  v'(x)  dx 


(H) 


Equation  11  is  the  form  actually  used  to  determine  the  function  x(p)  in  the  machine 
calculations.  This  relationship  may  alternatively  be  expressed  in  terms  of  fluxes 


-Bo/Bd  -  cos  a) 


Fo(0,L)  -  -V 

IT 


v'(x)  dx 


v/VBd 


(12) 


■cos 


o)  " : 


where 

v  v(x)  “  s/~K  F(B0/x,  L) 

L  is  held  a  constant  during  the  integration. 

Equation  12  multiplied  by  2w  has  been  machine  -programmed  to  calculate  the 
distribution  of  flux  $g|_(a)  per  unit  cosine  of  the  pitch  angle.  Thus  the  machine 
evaluates 


Bj/B  sin^a 


’(x) 


^Bq/B  sin^a  -  x" 


dx 


(13) 


Equation  13  has  one  interesting  immediate  application.  If  L  is  held  constant,  the 
numerical  value  of  the  flux  per  steradlan  about  some  pitch  angle  a^  at  a  location 
Bj  is  the  same  as  the  flux  per  steradian  measured  about  some  other  angle  a  at  another 
location  B,  if  a  and  B are  related  to  q^,  and  Bj  by 


.  2 
sin  a 


2 

sin  aj 


>1 
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In  other  words,  <f>gL(a)  ■*  determined  completely  by  knowledge  of  either  <^g^(ir/2) 

$q  .(a)  for  all  Land  a.  Again  B0  ■  B0(L)  is  the  equatorial 
Bo** 


for  all  B  and  L,  or  oy 
magnetic  induction. 


Some  plots  of  ^(a)  for  selected  L  values  are  shown  here  (Figures  16  through 
20)«  The  graphs  were  °  constructed  by  drawing  a  smooth  curve  through  points  calcu¬ 
lated  at  intervals  of  10  degrees  in  a.  From  these  graphs  the  pitch-angle  distribution 
anywhere  on  the  selected  L  values  may  be  determined. 


The  graphs  are  shown  for  the  four  maps:  electrons,  the  Hess  P-1  protons,  high- 
energy  protons,  and  low-energy  protons.  For  each  map  the  graphs  are  arranged  in 
order  of  increasing  L.  A  few  interesting  properties  of  these  curves  may  be  noted. 
Many  of  them,  as  shown,  have  a  nearly  linear  rise  in  the  plot  as  a  function  of  a. 

At  the  larger  L  values,  the  curves  seem  to  have  a  flat  top  or  plateau  with  some 
"wiggles.  "  Thus  the  curves  shown  for  L  ■  5  are  nearly  isotropic  except  for  a  cone  of 
about  10  degrees  half-angle  about  the  B  line.  The  behavior  of  the  curves  near  the 
loss  cone  angle  is  the  most  uncertain  feature,  but  this  depends  on  the  "edges"  of  the 
belt.  Work  should  be  done  on  a  description  of  the  edges  of  the  belts,  but  that  prob¬ 
lem  is  altogether  different  from  the  descriptive  type  of  problems  treated  here. 
Presumably,  a  better  understanding  of  the  loss  mechanism  is  required.  In  these  edge 
regions,  the  pitch-angle  distribution  is  not  so  simply  related  to  the  density  because 
scattering  affects  the  orbits  of  a  significant  fraction  of  the  particles  passing  through 
the  edges. 


These  curves  could  probably  be  improved  by  some  smoothing.  Many  of  them 
have  wiggles,  possibly  caused  by  inaccuracies  of  the  maps.  Already,  however,  some 
smoothing  has  occurred  In  their  construction.  New  flux  maps,  perhaps  more  accurate 
than  the  original,  could  then  be  calculated  with  Equation  9.  This  Integral  is  the 
inverse  of  that  used  to  determine  the  pitch-angle  distribution  from  the  flux  map. 

Since  the  Integral  (Equation  9)  and  its  Inverse  are  similar  in  form,  with  some  small 
modifications  the  angular  distribution  routine  could  also  be  used  to  calculate  its 
own  Inverse.  There  is  still  some  controversy  (Reference  30)  over  the  accuracy  of 
the  numerical  methods  used  to  evaluate  the  integral  (11).  That  is,  the  flux  maps 
have  been  smoothed  before  applying  the  transformation.  However,  one  may  apply  a 
matrix  transformation  (Reference  31)  to  approximate  the  integral.  The  matrix  may  be 
defined  so  as  to  operate  only  on  measured  flux  points.  This  can  lead  to  strange 
’bumps"  in  the  angular  distribution,  and  apparently  some  significant  deviations  from 
the  distribution  shown  here.  Further  investigation  of  the  validity  of  several  numeri¬ 
cal  procedures  would  be  worthwhile  for  application  as  better  flux  maps  become 
available,  but  the  approach  presented  here  seems  best  for  now. 

Furthermore,  from  the  calculated  pitch-angle  distributions,  better  geometrical 
factors  and  efficiencies  could  be  computed  for  counters.  In  this  light,  data  already 
gathered  could  be  reinterpreted  and  future  experiments  could  be  designed  for  better 
accuracy. 
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3.2  UNTRAPPED  RADIATION 


Galactic  cosmic  rays  and  high-energy  particles  ejected  from  the  sun  during 
large  solar  flares  comprise  the  untrapped  radiation  environment  in  the  vicinity  of  the 
earth.  Experimental  results  indicate  that  both  types  of  radiation  are  approximately 
isotropic  outside  the  magnetosphere.  That  is,  the  maximum  amplitude  of  the  average 
galactic  cosmic  ray  anisotropy  is  4  x  10"3  times  the  average  cosmic  ray  flux 
(Reference  32).  Solar  flare  particles  usually  reach  the  earth  nearly  Isotropically, 
within  5  to  10  percent  (References  33  and  34).  However,  during  a  short  time  at  the 
beginning  of  some  events  high  energy  particles  of  Intensities  10  to  100  times  the 
maximum  isotropic  radiation  intensity  do  arrive  with  significant  anisotropies  (Ref¬ 
erence  35).  At  the  present  stage  of  development,  the  insensitivity  of  total  mission 
dose  to  such  anisotropy  does  not  warrant  the  extensive  computer  analysis  required  to 
determine  the  effects  of  the  anisotropy  on  the  flux. 

The  method  adopted  for  computing  doses  from  untrapped  radiation  is  to  specify 
effective  spectra,  discussed  below,  fora  set  of  L  shells  and  to  determine  the  total 
flux  that  reaches  the  vehicle  in  each  L  shell.  The  dose  is  then  computed  from  a  set 
of  dose  conversion  numbers  calculated  by  the  shielding  codes,  using  as  input  the 
spectra  for  the  different  L  shells.  The  input  spectra  can  Include  both  sources  to¬ 
gether  if  desired.  However,  since  the  galactic  cosmic  radiation  constitutes  a  small 
contribution  to  the  dose  for  orbiting  vehicles,  and  since  it  changes  very  slowly  in 
time,  its  effect  is  easier  to  compute  separately,  as  is  discussed  in  Section  3.2.2. 

3.  2. 1  Solar  Particle  Events 

In  constructing  the  solar-flare  proton  subroutine  several  simplifying  assumptions 
have  been  made  about  the  time  development  of  the  energy  spectrum  of  the  protons 
outside  the  magnetosphere.  The  available  experimental  data  indicate  that  solar  par¬ 
ticle  events  have  the  following  general  features.  A  few  minutes  after  the  optical 
maximum  of  the  solar  flare,  the  first  high-energy  particles  reach  the  earth.  The 
particle  intensity  increases  as  the  lower-energy  particles  arrive  until  a  maximum  is 
reached.  The  time  of  the  maximum  Intensity  Is  from  8  to  30  hours  after  the  beginning 
of  the  event.  After  the  maximum,  the  Intensity  decreases  slowly  to  zero.  The  total 
duration  of  the  event  ranges  from  2  to  5  days.  Except  for  the  first  few  minutes  of  the 
event,  the  velocity  distribution  outside  the  magnetosphere  Is  roughly  isotropic.  In 
most  events  the  particles  reaching  the  earth  are  predominantly  protons,  although.  In 
a  few  events,  appreciable  numbers  of  alpha  particles  have  been  observed.  For  this 
subroutine,  it  will  be  assumed  that  all  the  particles  are  protons.  At  any  time  during 
the  event,  the  particle  spectrum  can  be  adequately  described  by  an  expression  of  the 
form 


J(>r)  -  i0w  .xpi-f/p0(t» 
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where  P  is  the  particle  rigidity  (momenturVchange),  j0(t),  and  PQ(t)  are  functions  of 
the  time  only,  and  j(>P)  Is  the  flux  of  particles  with  rigidity  >P.  The  rigidity  is  re¬ 
lated  to  the  kinetic  energy  E  of  the  particle  as  followst 

P(E)  -  2  J  E2  +  2  me2  E 

A  convenient  set  of  units  gives  P  In  Mv,  If  2  is  the  .number  of  electronic  charges,  me7 
is  the  particle  rest  mass  in  energy  units  (Mev),  and  E  Is  expressed  in  Mev. 

The  functions  jQ(t)  and  P0(t)  are,  in  general,  quite  complicated.  However,  it 
has  been  observed  that  the  integral  spectrum  for  the  whole  event  can  usually  be  ex¬ 
pressed  in  the  same  general  form.  Thus 


J(>P)  -  JQ  exp(-P/Pc)  -  j  jQ(t)  exp[-P/PQ(t)]  dt 


where  P0  is  roughly  equal  to  the  value  PQ(t)  at  the  time  of  maximum  flux.  To  simplify 
the  input  data  requirements,  the  following  assumptions  are  made  concerning  the  time 
behavior  of  jQ(t)  and  PQ(t). 


|  (t)  exp[-P/P Q(t)]  -  J0exp(-P/PQ)  F(t) 


where 


0  »  '  <  w  -  7tr 


1.001 


t^eXpC('*W)/TR'  '  'max  -  7t„  <  *  <  », 


'R  D 


max 


F(t) 


001 


Tr+t^  exp^max”^/TD^'  ^nxjx  -  *  -  ^x  +  7td 


0  ,  t  +  7Tn  <  t 
'  max  D 


and  t^,  Tq,  and  t  are  the  rise  time,  decay  time  and  time  of  maximum  intensity 
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for  the  event.  For  the  purpose  of  this  subroutine,  the  exponential  rise  and  decay  in 
intensity  will  adequately  predict  the  total  dose  and  dose  rate. 

The  effect  of  the  earth's  magnetic  field  on  the  Incident  radiation  must  be  in¬ 
cluded.  The  theory  of  the  allowed  cone,  as  developed  by  St&Vmer  (Reference  36), 
will  be  used  as  the  basis  for  calculating  the  modification  of  the  energy  spectrum. 

This  theory  is  based  on  a  dipole  approximation  to  the  geomagnetic  field.  The  con¬ 
stants  of  the  motion  of  a  charged  particle  in  a  dipole  field  are  the  energy  and  the 
canonical  momentum  about  the  dipole  axis.  With  these  constants,  the  differential 
equations  governing  the  trajectory  can  be  reduced  to  one  equation  of  motion  in  the 
meridian  plane  (References  37  and  38).  Stormer  obtains  an  expression  for  the  angle 
that  the  trajectory  makes  with  the  meridian  plane  and  deduces  that  the  trajectory 
must  lie  within  a  cone  of  allowed  directions  around  the  normal  to  the  meridian  plane. 
The  angle  of  the  cone  varies  with  the  energy  of  the  particle  and  with  the  position  of 
the  observation  point.  Lemaitre  used  Liouville's  theorem  to  prove  that  If  the  flux  of 
particles  is  isotropic  outside  the  magnetosphere,  then  the  intensity  in  a  given  direc¬ 
tion  within  the  allowed  cone  is  the  same  as  the  intensity  outside  the  magnetosphere. 
Hence,  the  energy  spectrum  inside  the  magnetosphere  is  given  by  that  outside  times 
the  ratio  as  a  function  of  energy  of  the  solid  angle  of  the  allowed  cone  to  4ir. 
Stormer  finds  that  the  allowed  cone  is  defined  by  the  equation 
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where  M  is  the  dipole  moment  of  the  earth  (0.311653  gauss  earth  radii  )  X  is  the 
geomagnetic  latitude,  R  is  the  Invariant  radius  in  earth  radii,  P  is  the  particle  rigi¬ 
dity  defined  above,  and  y  is  the  half  angle  of  the  allowed  cone  about  the  normal  to 
the  meridian  plane.  Thus,  we  see  that  the  cone  is  completely  closed  for  particles  of 
rigidity  leu  than 


P 


1.497x10 
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These  particles  cannot  reach  the  point  (R,  L)  from  outside  the  magnetosphere.  The  cone 
is  completely  open  for  particles  of  rigidity  greater  than 

p  .  1.497x104  (  2 

l5  \i  +  /i  -aw57* 

Figure  21  is  a  plot  of  these  two  expressions  at  an  altitude  of  300  nautical  miles  as  a 
function  of  L  From  the  expreuions  above,  it  is  easy  to  see  that  the  spectral 
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FIGURE  21,  Proton  cutoff  energies  versus  L 


52 


modification  is  most  strongly  dependent  upon  L  Thus,  we  assume  that  a  resultant 
spectrum  can  be  defined  for  any  L  shell  such  that  at  any  point  on  the  shell,  the 
instantaneous  spectrum  is  approximated  by 


J(  >P)  -  J0Gi(P,P0)  F(t) 

th 

where  the  Gj  are  the  modified  spectra  on  the  i  L-shell  and  F(t)  describes  the  flux¬ 
time  dependence.  The  flux  for  a  specific  mission  is  calculated  as  follows.  The  L 
range  is  divided  into  7  macro-shells  such  that  the  spectral  modification  is  relatively 
slight  within  each  macro-shell.  The  values  L.  will  be  used  to  designate  the  interval 
L.  i  <  L  <  L. .  The  seven  L's  are  given  as  Input  data.  Since  the  trajectory  code 
calculates  L  as  a  function  of  time,  it  is  possible  to  obtain  a  time  integration  of  the 
function  F(t)  for  that  portion  of  the  mission  spent  in  each  of  the  L  intervals.  The 
total  flux  reaching  the  vehicle  is  then 
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where  Fj  is  the  time  integral  of  F(t)  over  the  time  spent  in  the  corresponding  L  interval 

In  addition  to  calculating  the  time  integral  of  F(t)  the  subroutine  also  keeps  a 
running  total  of  the  dose  received  in  the  i™  L  interval  by  using  a  set  of  dose  conver¬ 
sion  numbers  k;  provided  by  the  shielding  code. 

The  subroutine  has  been  designed  for  modified  exponential  rigidity  spectra. 

Since  the  functions  Gj  must  be  specified  as  input  data,  other  functional  forms  may  be 
used  such  as  power  laws  in  energy.  The  suggested  modification  of  the  exponential 
rigidity  spectrum  is  an  abrupt  cutoff  below  the  vertical  cutoff  rigidity.  Thus, 

fo  ,  P  <  PC(L,) 

G|(PrP„>  -{ 

/  «xp(-p/Po),  p  >  pea,) 

where 

Pc(l,)  -  1.4?7(1oVl,2 

The  geomagnetic  storm  that  accompanies  the  larger  events  produces  a  reduction 
in  the  strength  of  the  magnetic  field  so  that  the  effective  L  value  at  a  given  point  Is 
larger  than  the  normal  L  value  (see  Figure  22).  This  modification  of  the  geomagnetic 
field  can  be  Included  in  the  spectra  Gj  by  using  the  effective  L  rather  than  the  normal 
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Proton  cutoff  energy  (Mev) 


L  in  the  definition  of  Pc.  For  example,  the  set  of  Lj  given  by  *  3.0,  Lj  *  3.6, 

L3  ■  4. 2,  L4  *  4. 7,  L5  *  5. 2,  =  5.6,  and  L7  =  00  would  correspond  to  proton 

vertical -cutoff  energies  in  Mev,  of  E^  "610,  E2  *  241,  E3  ■  106,  £4  =  58. 3, 

E5  ■  32.5,  E^  *  20.5,  and  E7  ■  0  during  the  event.  During  quiet  times  the  corres¬ 
ponding  vertical  -cutoff  energies  would  be  E]  ■  980,  E2  =  560,  E3  ■  325,  E4  *  217, 

E5  -  150,  E6  =»  114,  E7«0. 

The  accuracy  of  the  experimental  data  on  fluxes  and  energy  spectra  does  not 
justify  a  completely  rigorous  treatment  of  the  problem.  The  L-ehell  approximation 
could  be  improved  upon  by  a  method  that  computes  the  modified  spectrum  at  each 
point  on  the  trajectory.  It  is  recognized  that  the  L-shell  approximation  is  a  first 
attempt  to  compute  the  effect  of  the  geomagnetic  field  on  untrapped  radiation  and 
does  have  the  advantage  of  fitting  easily  into  the  calculational  procedures  established 
for  the  trapped  radiation. 

3.2.2  Galactic  Cosmic  Rays 

The  radiation  dose  to  space  vehicles  from  the  galactic  cosmic  radiation  can  be 
calculated  from  the  measurements  of  Neher  and  Anderson  (References  39  and  40)  and 
Lin,  et  al.  (Reference  41).  Neher  and  Anderson  have  flown  ion  chambers  in  balloons 
at  polar  latitudes  for  a  number  of  years.  They  flew  a  similar  ion  chamber  aboard  the 
Mariner  2  space  probe,  and  a  comparison  of  the  Mariner  2  and  the  balloon  measure¬ 
ments  (Reference  40)  indicated  that  21  percent  of  the  dose  to  the  balloon-borne 
instruments  was  due  to  albedo  arising  from  the  interaction  of  the  primary  radiation  with 
the  earth's  atmosphere  and  magnetic  field.  These  measurements  can  be  expressed  as 
dose  in  rad  hr  (instead  of  the  ions  cm”3  iec~*  Atm”*  of  air  as  reported  by  the 
authors).  The  values  used  in  the  present  treatment  include  the  albedo  contribution, 
and  are  for  particles  incident  from  the  upper  hemisphere;  this  treatment  is  appropriate 
for  vehicles  in  near-earth  orbits,  but  for  vehicles  a  large  distance  from  the  earth 
these  values  should  be  reduced  by  20  percent,  and  then  corrected  for  spherical  geo¬ 
metry.  Figure  23  shows  the  variation  of  ion  chamber  response  with  atmospheric  depth 
near  the  magnetic  pole  for  several  times  during  the  solar  cycle,  with  a  scale  added  to 
indicate  the  dose  rate  in  mrac|/hr.  It  can  be  noted  that  the  dose  changes  very  little 
with  atmospheric  depth  near  the  top  of  the  atmosphere;  thus  the  maximum  dose  values 
can  be  used  to  evaluate  the  dose  inside  space  vehicles  for  shielding  thicknesses  cur¬ 
rently  under  consideration.  It  can  also  be  seen  that  the  dose  rate  changes  by  a  factor 
of  two  during  the  sunspot  cycle,  the  maximum  dose  rate  from  galactic  cosmic  radiation 
coming  at  the  same  time  as  the  minimum  in  solar  activity. 

These  measurements  of  the  dose  rate  can  be  related  to  the  measurements  of  Lin, 
et  al.  (Figure  24).  These  results  indicate  that  the  counting  rate  from  galactic  cosmic 
rays,  as  measured  by  the  Explorer  VII  satellite,  is  constant  above  L  ■  2.5  and  is 
approximately  a  linear  function  of  L  between  L  values  of  1.1  and  2.5.  Hence,  the 
dose  rate  will  also  depend  linearly  on  L,  provided  that  the  dose  per  particle  is 
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FIGURE  24.  Cosmic 'ray  counting  rate  as  measured  by  Explorer  VII  satellite  (November  1959  to  December  1960) 


constant.  However,  the  variation  in  dose  per  particle  is  such  that  estimates  based  on 
this  approximation  are  conservative:  the  dose  is  slightly  overestimated  in  equatorial 
regions,  where  the  dose  is  small  anyway  due  to  a  smaller  particle  flux. 

Normalizing  the  dose  rate  measurements  of  Neher  and  Anderson  to  the  L  varia¬ 
tion,  we  have  the  dose  rate  D  as  a  function  of  L 


r[1. 92(L-1)  +  0. 72] 
6(10  7)  miad/hr. 


(10  7)  mrad/hr. 


for  1  <  L  <  2.5 
forL  >2.5 


at  the  time  of  maximum  in  the  solar  cycle. 


For  the  time  of  solar  minimum 


C4.32(L-1)  +  0.72]  (10-7)  mrad/h 
7.2(10“7)  m rad/hr. 


for  1  <L  <2.5 
forL  >2.5 


The  resulting  dose  rates  will  be  much  smaller  than  those  arising  from  the  trapped  par¬ 
ticles  for  most  cases  of  interest.  The  maximum  galactic  dose  rate  is  17  millirad/day 
during  solar  minimum  and  9  millirad/day  at  the  time  of  solar  maximum.  For  cases 
where  the  dose  from  trapped  particles  exceeds  these  values  significantly,  it  is 
recommended  that  the  values  for  L  >2. 5  be  used  to  estimate  the  galactic  radiation 
dose. 
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4.0  VEHICLE  ANALYSIS 


An  evaluation  of  the  shielding  of  a  space  vehicle  for  any  Incident  radiation 
flux  Is  dependent  on  an  analysts  of  the  geometry  and  material  composition  of  the 
vehicle.  This  analysts  Is  designated  as  a  "sector  analysts"  of  the  vehicle,  and  Is 
accomplished  by  placing  the  origin  of  a  spherical  coordinate  system  at  the  dose  point 
and  determining  the  shield  thickness  of  the  vehicle  as  a  function  of  solid  angle.  The 
vehicle  Is  divided  Into  Nsectors,  each  subtending  a  solid  angle  On  where 

N 

Z  "n  ‘  4" 

n*1 

The  slant  thickness  of  the  materials  providing  shielding  In  a  sector  can  be  determined 
as  Figure  25  will  show.  Mathematical  refinements  can  be  made  to  improve  the 
accuracy  of  the  analysts  by  analyzing  many  vectors  per  sector  and  determining  some 
effective  average  thickness.  Boeing  and  other  organizations  (References  42  and  43) 
have  been  developing  computer  programs  to  perform  such  a  geometric  analysis  auto¬ 
matically.  The  self-shielding  of  the  human  body  can  also  be  Included  In  the  analysis 
when  the  dose  point  Is  In  an  astronaut. 

In  the  original  Boeing  space-radiation  programs  the  dose  penetrating  the  vehicle 
was  computed  by  the  expression 

N 

D,(p)  *  S!  TIn(P}  wn 
n"I 

where  D.(p)  Is  the  dose  at  a  point  p  from  radiation  of  type  I;  Sj  Is  the  flux  of  radiation 
of  type  I;  Tjn(p)  Is  the  dose  transmission  factor  for  radiation  type  I  penetrating  sec¬ 
tion  n,  and  un  is  a  weighting  factdr  defined  as  0^/4 ir,  where  Qn  Is  the  solid  angle 
subtended  by  section  n. 

The  data  necessary  to  compute  Tjn(p)  is  the  energy  spectrum  of  the  incident 
radiation  and  the  atomic  number  and  thickness  of  each  material  in  a  given  sector 
(ordered  from  Incident  face  to  exit  face). 

To  improve  the  mathematical  model  to  include  the  effects  of  the  angular  distri¬ 
bution  of  the  penetrating  radiation  and  the  scattering  within  a  vehicle,  more  Input 
data  would  be  required*  A  simplified  sector  of  a  vehicle  and  the  variables  that  are 
required  for  an  Improved  radiation  analysis  are  shown  in  Figure  25.  To  achieve  an 
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Required  data  for  each  layer,  if  in  sector  n 
Z.^  *  atomic  number  of  layer  (I,  n) 

Rjn  *  scalar  distance  between  dose  point  and  near  surface  of  layer  (ir  n) 

R  jo  *  normal  to  layer  (i, n) 

t.^  *  thickness  of  layer  (I,  n) 

p  »  demfty  of  layer  (i,  n) 

6,  *  angle  between  N,  and  R, 
in  In  in 

*  solid  angle  subtended  by  sector  n 

?  “  vector  of  axis  of  azimuthal  symmetry  for  Incident  flux 

-*  — * 

a  m  angle  between  R,  and  I 
n  *  In 

FIGURE  25.  Input  data  for  improved  vehicle  analysis 
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improved  analysis,  eight  variables  for  each  layer  in  each  sector  must  be  specified 
compared  to  the  three  variables  required  for  the  original  analysis.  The  following 
sections  of  this  report  describe  the  improved  computations  of  the  dose  transmission 
factor  using  either  3  or  8  input  variables  for  various  types  of  incident  radiation  and 
their  secondaries.  It  will  be  established  that  the  use  of  an  improved  analysis  em¬ 
ploying  the  eight  variables  is  not  justified  at  this  time  because  the  small  increase  in 
accuracy  does  not  warrant  the  extra  effort  required  to  obtain  the  input  data  and 
perform  the  additional  analyses. 
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5.0  ELECTRON  SHIELDING 


The  prediction  of  dosage  within  a  space  vehicle  exposed  to  electron  radiation 
is  a  complex  problem  influenced  by  many  variables.  The  purpose  of  this  study  was  to 
examine  variables  such  as: 

1)  The  incident  angular  and  energy  distributions; 

2)  The  LET  and  energy  spectra; 

3)  The  angular  distribution  of  the  penetrating  radiation; 

4)  The  use  of  multislabs  of  various  materials;  and 

5)  The  details  of  input  geometric  sectoring  for  the  vehicle  and  crew. 

The  prior  Boeing  electron  shielding  codes  (References  44,  45,  and  46)  were  used 
In  some  of  the  analyses  and  formed  the  basis  for  the  codes  developed  in  this  study.  In 
many  situations,  these  codes  are  satisfactory  for  space  shielding  calculations. 

5.1  MONTE  CARLO  STUDIES 

5.1.1  Interactions 

The  object  of  this  study  was  to  develop  a  mathematical  model  to  compute  the 
transmitted  number,  energy  spectrum,  and  angular  distribution  of  electrons  penetrating 
a  shield.  An  electron  penetrating  a  shield  will  follow  a  tortuous  path  as  it  interacts 
with  atomic  electrons  and  nuclei  of  the  shield.  Even  though  the  total  path  length  of 
monoenergetic  electrons  Is  relatively  constant,  the  distance  between  the  face  of  the 
shield  and  the  end  of  an  electron  path  will  vary  greatly  due  to  scattering.  This  fact 
makes  the  electron  range  a  poor  parameter  in  electron  shielding  calculations. 

Electrons  In  a  shield  can  be  scattered  by  four  principal  interactions: 

1)  An  inelastic  collision  with  atomic  electrons,  where  one  or  more  atomic 
electrons  are  raised  to  an  excited  state  or  unbound  state.  This  Is  the  predominant 
scattering  mechanism  for  electrons  with  a  few  hundred  kev  kinetic  energy; 

2)  An  inelastic  collision  with  a  nucleus,  where  the  electron  is  deflected  and  a 
quantum  of  radiation  is  emitted  (bremsstrahlung); 

3)  An  elastic  collision  with  a  nucleus,  where  the  electron  loses  only  the  kinetic 
energy  necessary  to  consorve  momentum  between  the  two  particles;  and 


4)  An  elastic  collision  with  the  atom,  only  of  Importance  for  electrons  with 
energies  less  than  160  ev. 

Detailed  descriptions  of  these  electron  interactions  are  presented  by  Evans 
(Reference  47)  and  Birkoff  (Reference  48). 

The  Boeing  electron  Monte  Carlo  code  (BEMC)  (Reference  45)  was  developed 
using  the  computational  models  of  Perkins  (Reference  49)  and  Leiss,  Penner,  and 
Robinson  (Reference  50).  The  complicated  process  of  electron  interaction  in  a  shield 
is  simplified  into  two  mechanisms:  radiation  loss  and  ionization  loss.  An  electron  is 
followed  through  a  series  of  thin  slabs  where  it  is  assumed  to  have  suffered  discrete 
energy  losses  and  angular  deflection.  In  each  increment  the  energy  loss  due  to 
Ionization  was  obtained  from  the  average  ionization  loss  tables  of  Nelms  (Reference 
51).  As  shown  in  Figure  26,  the  ionization  energy  loss  in  aluminum  changes  rapidly 
with  energy  below  1  Mev.  Therefore,  the  slab  increments  should  decrease  with 
decreasing  electron  energy  below  1  Mev  so  that  the  loss  can  be  described  by 

Elon  lost  * 


where  Ax  is  the  increment  thickness  and  ^(E)  i*  the  energy  loss  for  an  electron  of 
energy  E.  The  machine  program  permits  the  use  of  three  forms  of  Ax:  a  constant 
value,  a  constant  multiplied  by  the  electron  energy,  and  a  constant  divided  by  the 
electron  energy.  Values  of  Ax  from  0.025  to  0. 1  of  the  electron  range  were  tested 
and  found  to  provide  answers  within  the  statistical  error  associated  with  the  number  of 
histories  processed.  Minimum  deviations  were  observed  using  Ax/E,  where  Ax  *  0.05  R 
(R  is  the  range  in  grr/cm^  and  E  is  In  Mev). 


The  second  criterion  for  the  increment  thickness  was  that  a  small  net  angular 
deflection  of  the  electron  be  obtained  In  the  increment.  For  the  BEMC,  the  scatter¬ 
ing  angle  is  determined  by  the  Bethe  modifications  to  the  Moliere  relationships  (Ref¬ 
erence  52)  which  are  only  valid  for  small  scattering  angles  (less  than  20  degrees).  A 
thickness  of  1/20  of  the  electron  range  Is  normally  used  for  Ax.  (See  Reference  44  for 
full  details).  A  study  of  the  current  multiple -scattering  theories  for  electron  penetra¬ 
tion  was  conducted  and  the  analysis  Is  presented  In  the  Appendix.  The  Goudschmldt- 
Saunderson  (Reference  53)  multiscattering  theory  was  found  to  be  superior  at  large 
angles;  however,  the  computation  time  required  on  a  computer  Is  excessive.  This 
method  has  been  employed  by  Berger  (Reference  54)  in  an  electron  Monte  Carlo  study 
and  reasonable  agreement  is  obtained  between  his  data  and  the  BEMC  data  using  thin 
slabs. 

The  probability  of  radiation  loss  P(E-£)  Is  obtained  from  the  relationships 
derived  by  Bethe  and  Heitler  (Reference  55)  expressed  as 
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FIGURE  26.  Electron  energy  los*  in  aluminum  (Reference  51) 
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ln(E0/E)[(t/|n2)-l)3 

T  Ct/ln  2) 


P(E0-E)  = 


In  this  equation  E0  is  the  initial  electron  energy,  E  is  the  final  energy  after  radiation 
loss,  and  t  is  the  ratio  of  the  slab  thickness  to  the  radiation  length  XQ  where 

_L  *  z2  r02  in(183Z“1/3) 

-13 

where  a  *  1/137,  N0  is  Avogadro's  number,  and  r0  M  2. 28  x  10  cm 

The  Monte  Carlo  method  neglects  straggling  and  does  not  follow  delta  or 
knock-on  electrons.  Comparison  of  BEMC  data  with  existing  information  shows  that 
these  phenomena  are  not  important  to  the  primary  electron  number  density  or  angular 
distribution.  These  corrections  were  made  by  Perkins  (Reference  49)  in  his  Monte 
Carlo  study  and  were  found  to  be  negligible. 

5.1.2  Modifications 

The  original  BEMC  was  written  In  machine  language  and  was  converted  to  For¬ 
tran  for  compatibility  with  the  AFWL  computers.  Modifications  were  made  to  provide 
the  computation  of: 

1)  The  cosine  of  the  exit  angle  of  the  electron; 

2)  The  summation  of  number  spectra  over  angle; 

3)  The  summation  of  number  spectra  over  energy; 

4)  The  energy  spectrum  of  the  bremsstrahlung;  and 

5)  The  LET  number  spectrum. 

The  modified  code  will  accept  (1)  any  spectrum  of  electron  energies  including 
a  monoenergetic  spectrum  treated  as  a  special  case,  (2)  a  monodirectional  or  isotropic 
distribution  in  incident  angle,  and  (3)  a  combination  of  twelve  shielding  materials. 

The  code  will  compute  either  number  and  energy  spectra  for  various  exit  angles  or 
LET  spectrum  for  various  exit  angles,  and/or  bremsstrahlung  production.  The  changes 
are  described  in  the  operator's  manual.  Volume  II  of  this  report. 

5.1.3  Angular  Distribution  of  Penetrating  Electrons 

Modifications  were  made  in  the  treatment  of  the  angular  distribution  to  examine 
the  scattering  of  electrons  penetrating  a  shield,  to  obtain  a  model  for  scattering  inside 
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of  a  vehicle,  and  to  establish  the  angular  distribution  of  the  penetrating  electrons. 
The  average  cosines  were  calculated  for  each  increment  of  a  target  slab  according 
to  the  formula 


<  cos  0>  *  )>  nj  cos  6-/)>  nj 

•  • 

i  i 

where  n«  is  the  number  of  electrons  scattered  out  of  the  increment  at  the  angle  0j. 

The  angular  distributions  in  the  forward  direction  were  computed  in  increments 
15  degrees  wide  for  an  initial  flux  of  103  electrons  incident  at  0,  30,  and  60  degrees 
with  energies  of  1,  3,  and  6  Mev.  In  order  to  cover  a  wide  range  of  Z  uniformly, 
aluminum,  iron,  silver,  and  lead  were  used  as  target  materials. 

A  general  analysis  of  the  Monte  Carlo  electron  angular  distributions  is  not 
attempted  here,  since  no  completely  rigorous  multiple-scattering  theory  exists  which 
takes  into  account  all  the  effects  of  energy  loss,  straggling,  and  various  boundary 
conditions.  The  Monte  Carlo  method  itself  is  the  only  technique  available  at  present 
for  treating  a  general  class  of  problems. 

The  following  analysis  is  based  on  a  consideration  of  the  angular  distribution  for 
target  thicknesses  comparable  to  the  range  of  the  incident  electrons.  It  has  been 
previously  hypothesized  that  the  electron  angular  distribution  becomes  essentially 
isotropic  after  passage  through  sufficiently  thick  targets.  This  assertion  has  been 
verified  theoretically  by  Rohrlich  and  Carlson  (Reference  56)  for  the  simplest  case  of 
an  infinite  medium.  Obviously,  if  this  also  held  true  for  bounded  media,  the  prob¬ 
lem  of  angular  distributions  could  be  effectively  eliminated  by  taking  sufficiently 
thick  targets.  However,  for  plane  parallel  targets  of  thickness  comparable  to  the 
total  range,  the  experimental  results  obtained  by  Frank  (Reference  57)  indicate  that 
the  angular  distribution  is  not  isotropic,  but  remains  sharply  peaked  in  the  forward 
(and  backward)  direction.  In  addition,  Frank  states  that  after  the  electrons  have 
reached  a  certain  "diffusion"  depth  aQ,  the  half-width  of  the  angular  distribution 
remains  constant  so  that  further  penetration  by  the  electrons  may  be  dependent  only 
on  the  transmission  factor.  In  accordance  with  these  results,  the  Mon te  Carlo  data 
was  first  analyzed  to  determine  the  diffusion  lengths  and  then  the  actual  angular  dis¬ 
tributions  were  compared  at  cr0. 

The  diffusion  depth  a  for  each  energy  and  target  material  was  determined  from 
the  Monte  Carlo  data  by  estimating  the  point  at  which  the  average  cosine  of  the 
scattering  angles  for  normal  incidence  appeared  to  reach  its  asymptotic  value.  The 
"asymptotic"  point  was  assumed  to  be  equal  to  a  .  Ibis  assumption  was  justified  by 
the  agreement  in  the  asymptotic  <cos6>  values  for  Frank's  results  and  the  Monte 
Carlo  data.  The  average  cosine  is  defined  by 
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<co»  9>  ■  jco*  ®  9)  dQ/  ^f(0/9)  cl 0 

Frank's  experimental  angular  distributions  are  accurately  described  at  e0  by  the 
distribution  function  f(9)  ■  f(0)  cos ^9,  regardleu  of  the  target  material.  This  gives 
an  average  cosine  of  0.75.  The  Monte  Carlo  results  for  all  materials  and  energies 
Investigated  Indicated  a  value  of  <cos  8>  ■  0.74  ±  0.03  which  agrees  within 
statistical  error  with  Frank's  data  —  If  the  "asymptotic"  point  ts  Identified  with  the 
diffusion  length  aQ.  It  should  be  noted  here  that  the  assumption  of  an  Isotropic 
distribution  gives  <cos  9  >  ■  J^3. 

The  average  cosines  for  aluminum  at  1  Mev  for  Incident  angles  of  0,  30,  and 
60  degrees  are  shown  In  Figure  27.  The  error  brackets  on  the  curve  fcr  normal  Inci¬ 
dence  were  determined  by  the  statistical  error  In  the  number  of  electrons  transmitted 
through  each  Increment.  The  diffusion  lengths  for  aluminum.  Iron,  silver,  and  lead 
at  1,  3,  and  6  Mev  are  shown,  with  estimated  errors  In  Table  3. 


Table  3 

MONTE  CARLO  DIFFUSION  LENGTHS 
(c0  In  g/cm2) 


Element 

(Z) 

Initial  Energy 

1  Mev 

3  Mev 

6  Mev 

Al 

(13) 

0.17  ±  0.02 

0.85  *  0.10 

1.90  *  0.10 

Fe 

(26) 

O.il  ±  0.01 

0.55  *  0.08 

1.65  ±  0.10 

Ag 

(47) 

0.08  ±  0.02 

0.48  ±  0.05 

1.35  *  0.10 

Pb 

(82) 

0.075*0.01 

0.40  *  0.05 

0.95  *  0.10 

Since  Frank's  dataware  taken  only  at  1.75  Mev,  the  diffusion  lengths  can  be 
approximately  compared  through  the  modified  Moltere  theory  which  Frank  used  to  ob¬ 
tain  c0.  Frank  defined  the  condition  of  "total  diffusion"  to  occur  at  the  point  where 
the  angular  half-width  0h  was  equal  to  45  degrees.  A  good  approximation  for  0^, 
according  to  Hanson,  Lanzl,  Lyman  and  Scott  (Reference  58),  Is 


0H  -  0.833  XCV  B  -1.2  (14) 

which  Is  obtained  by  considering  only  the  dominant  Gaussian  term  of  the  angular  distri¬ 
bution.  According  to  Frank,  this  formula  holds  for  half-widths  In  excess  of  20  degrees. 
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FIGURE  27.  Average  cosines  for  1-Mev  electrons  incident  at 
0,  30,  and  60  degrees  on  aluminum 


The  quantities  Xc  and  B  are  given  in  the  Moliere-Bethe  theory  (Reference  52)  by 


Xc2  «  4irNde4  Z(Z  +  l)/(pv)2 
B  -  In  B  -  b 


05) 


where 

,  _  ,  r  6680pd  (Z  +  DZ1/3 
b  *  In  L  *  fy 

A(1  +  3. 34  aj) 


ai  =  ze2/*v 
N  *  atom  number  density 
A  *  atomic  weight 
v  *  electron  velocity 
B  ■  v/c 
P  *  momentum 


p  *  target  density 


These  equations  are  coupled  In  the  thickness  d.  To  find  the  diffusion  length  oq  s  pd 
for  8j»j  "45  degrees,  we  combine  Equations  Hand  15. 


Then 


9.58  8H2  A(py)2 


>H  'Hpv/ 


. « 90  r  a  4T  Cm)  1 

5,90  [  Z(Z  +  1)  (B -1.2)  1+4T  (T+l)] 

2 

where  T  is  the  kinetic  energy  in  Mev  of  the  incident  electrons  and  mQc  is  set  equal 
to  0.5.  Designating 


70 


f  -  5  90  [  _A  _4T2(T-H)  . 

f  5,90  L  Z{i  +  ])  1  +  4T(T+1)  ] 

(Z  +  1)Z1/3 

9  “  - 2" 

A(1  +  3. 34  ay 

the  resulting  expressions  for  o0  become 
oQ  -  f/(B  -  1.2) 


(16) 


B-lnB=ln[  ]  (17) 

r(B-1.2) 

For  specified  values  of  Z  and  T,  Equation  17  is  solved  for  B  by  iteration  and  aQ  is 
then  obtained  from  Equation  16.  The  calculated  values  of  a0  can  be  expected  to 
agree  only  qualitatively  with  the  Monte  Carlo  values  since  the  effects  of  energy  loss 
have  been  neglected.  The  results  are  shown  in  Figures  26,  29,  and  30  for  the  ener¬ 
gies  1, 3,  and  6  Mev,  respectively.  The  values  given  by  Frank  at  1.75  Mev  are  also 
shown  in  Figure  2B.  The  qualitative  agreement  is  good,  although  large  deviations 
occur  for  high  and  low  Z  values. 

The  diffusion  length  is  measured  along  the  incident  electron  direction  and  is  ob¬ 
served  to  increase  slightly  with  angle  measured  from  the  normal.  The  diffusion  lengths 
for  the  30°  and  60°  curves  shown  in  Figure  27  are  obtained  by  dividing  the  normal 
areal  density  by  the  cosine  of  the  angle. 

The  angular  distributions  of  transmitted  electrons  were  tabulated  as  histograms  with 
15-degree  increments.  Figure  31  shows  an  example  for  1-Mev  electrons  incident  on 
aluminum  of  thickness  0»175gm/cir  .  The  distribution  appears  to  be  best  described 
by  Frank's  empirical  cosine -squared  law.  The  function  (0.717  cos 0  +  cos  20)  sin  6 
obtained  by  Bethe,  Rose,  and  Smith  (Reference  59)  neglecting  energy  loss  effects.  In 
a  previous  Monte  Carlo  calculation,  Berger  (Reference  54)  obtained  a  better  fit  to 
his  data  with  this  function  than  with  the  cosine -squared  law.  However,  Berger's  dis¬ 
tribution  was  calculated  for  a  lower  energy  of  0.5  Mev  which  apparently  resulted  in 
an  increase  of  large-angle  scatterings. 

The  value  of  the  average  cosine  at  oQ  also  indicates  that  the  cosine -squared  law 
provides  a  better  description  of  the  angular  distribution.  The  function  (0.717  cos0  + 
cos2®)  sin0  gives  <cos  0>  ■  0.706  compared  with  0.750  for  Frank's  cosine-squared 
law  and  0.74  ±  0.03  for  the  Monte  Carlo  calculation.  For  target  thicknesses  greater 
than  a0,  the  cosine -squared  law  also  provides  a  reasonable  fit  as  shown  in  Figure  32. 
For  pd>  2o0,  the  results  are  uncertain  due  to  the  large  statistical  errors  which 
occurred  in  the  Monte  Carlo  data.  For  the  case  shown  in  Figure  30,  the  angular 
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FIGURE  28. 


Diffusion  lengths  for  1-Mev  electrons,  normal  Incidence 
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FIGURE  29.  Diffusion  lengths  for  3-Mev  electrons,  normal  incidence 
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FIGURE  31. 


Angular  distribution  of  1-Mev  electrons 
transmitted  through  0. 175  gm/cm^  of 
aluminum,  normal  incidence 
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FIGURE  32.  Angular  distribution  of  1-Mev  electrons 

transmitted  through  aluminum,  normal  incidence 
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increment*  at  2.2  a0  each  contained  less  than  40  electron  histories,  resulting  in  a 
statistical  error  greater  than  15  percent. 

From  the  preceding  considerations,  the  Monte  Carlo  results  indicate  that  the 
electron  angular  distributions  can  be  approximately  described  by 

f(0)  dfi  *  fQ  cos^9  sin9  d9  d$  (18) 

in  the  range  aQ  <  dp  <  2a0 .  The  normalizing  factor  fQ  is  a  function  of  the  thick¬ 
ness,  the  incident  energy  and  flux,  and  the  target  material.  It  can  be  determined 
from  the  Monte  Carlo  data  by  fitting  Equation  18  at  the  peak  angle.  The  results  for 
1-Mev  normal  incidence  with  an  initial  flux  of  10^  electrons  are  shown  in  Figure  33. 
For  all  materials  investigated,  f0  varies  linearly  with  thickness.  The  two  deviating 
points  of  the  iron  Monte  Carlo  calculations  are  considered  uncertain  due  to  the  rela¬ 
tively  poor  curve  fit  obtained.  However,  they  could  easily  have  been  made  col  inear 
with  the  other  points  if  the  distributions  within  the  angular  increments  used  were  ap¬ 
propriately  skewed.  The  linear  behavior  of  f0  is  associated  with  the  approximately 
linear  behavior  of  the  transmission  factor  in  the  range  of  thicknesses  considered.  The 
relation  between  fQ  and  the  transmission  factor  can  be  used  to  obtain  values  of  fQ  for 
different  energies,  incident  fluxes,  and  target  materials.  Since  f0  is  independent  of 
the  scattering  angle,  it  is  proportional  to  the  transmitted  flux  and  thus  to  the  trans¬ 
mission  factor.  Consequently,  the  ratio  of  different  values  off  is  equal  to  the  ratio 
of  the  appropriate  transmission  factors.  This  proportionality  may  also  justify  extending 
Equation  18  as  far  as  the  extrapolated  range  of  the  electrons. 

Analysis  of  the  Monte  Carlo  electron  angular  distributions  has  indicated  that  the 
Monte  Carlo  method  is  essentially  in  agreement  with  experimental  results.  The  hypo¬ 
thesis  of  isotropy  has  been  shown  to  be  invalid  for  thick,  bounded  media. 

The  behavior  of  the  angular  distribution  for  thicknesses  greater  than  oQ  can  be 
described  by  Equation  18.  The  general  transmission  factor  can  then  be  written  as  the 
product  of  the  distribution  function  f(0)  and  the  usual  angle-independent  factor  f  . 
Utilization  of  this  distribution  will  be  discussed  in  Section  7.0,  the  analysis  of  the 
importance  of  angular  distribution  of  incident  and  penetrating  elections. 

5. 1.4  Accuracy  of  Monte  Carlo  Calculation 

The  two  principal  sources  of  error  in  the  Monte  Carlo  treatment  of  electron  multi¬ 
ple  scattering  are  the  statistical  fluctuations  and  the  basic  assumptions  of  the  under¬ 
lying  theory.  The  statistical  error  can  be  made  arbitrarily  small  by  increasing  the 
sample  size;  however,  the  errors  caused  by  assumptions  and  approximations  in  the 
theory  are  often  neglected  and  may  lead  to  difficulty  if  the  theory  is  used  beyond  its 
range  of  validity.  This  is  particularly  true  in  multiple  scattering  where  the  require¬ 
ment  of  small  scattering  angles  is  usually  imposed.  To  estimate  the  magnitude  of  such 
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FIGURE  33.  Normalizing  factor,  f  ,  for  1000  electron 
histories  (normal  incicfenceat  1  Mev) 
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"theoretical  "  errors,  an  approximate  analysis  is  undertaken  here  of  the  consequences 
involved  in  applying  small -angle  theories  in  the  Monte  Carlo  method. 

The  Monte  Carlo  treatment  of  multiple  scattering  is  considered  valid  for  all 
scattering  angles  because  of  the  assumption  that  all  large-angle  scatterings  are  the 
result  of  many  successive  small -angle  scatterings.  Since  the  target  is  usually  divided 
into  small  increment  slabs  in  the  Monte  Carlo  method,  this  assumption  requires  that 
the  angle  of  scattering  from  each  slab  be  equal  to  or  less  than,  the  theoretical  limit. 

The  "theoretical"  error  can  thus  be  obtained  by  considering  each  slab  individually. 

For  increments  much  smaller  than  the  electron  range,  the  energy  loss  may  be 
neglected.  For  this  to  be  valid  for  all  target  increments,  the  Increments  may  have  to 
be  taken  successively  thinner  for  very  thick  targets. 

In  a  beam  of  electrons  pawing  through  a  single  increment,  nQ  electrons  will  make 
collisions  for  which  the  scattering  angle  9  exceeds  some  specified  maximum  scatter¬ 
ing  angle  •  If  it  is  assumed  that  any  single  collision  is  independent  of  the  others, 
then 

nQ  “  2tt  N  $  d  ^  a  (9)  sin9  d  9  (19) 

®max 

where  N  Is  the  number  of  target  atoms  per  cm 3,  d  is  the  Increment  thickness  and  <$> 
is  the  incident  flux.  For  the  screened  Coulomb  potential  V  *  (Ze/r)  exp(-r/a),  the 
scattering  cron  section  in  the  first  Bom  approximation  is  given  by 


<K0)  »  2  i ■■■  " - j 

p  v  (1  -cos9  +  26) 

where  6  ■  ti2/4a2p2,  and  p  and  v  are  the  momentum  and  velocity,  respectively,  of 
the  scattered  electron.  Performing  the  integration  In  Equation  19,  n0  Is  given  by 


Z2  4 

nQ  -  2w  N  *  d  — £-5 


+  28  -  cos9max)"1 


-[2(1  +8)]” 


Using  the  relation 

pV  “ 


4T2(T  +  l)2 
i+4t(t+iy 
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where  mQcz  **  1/2  and  T  is  the  kinetic  energy  in  Mev, 


,  /n\  k  i  .,-,2  2  r  1  +  4T(T  +  1), 
nQ  =  vr/8)  N  <1>  d  Z  r0  L - = - *-J 

T2(T+ir 


(1+28  -cose^J"1  -C2(l+6)]_1 


] 


2  . 


where  rQ  *  e  /mQc  is  the  classical  electron  radius.  The  quantity  8  is  determined 
from  the  screening  length 

,-1/3 


a  =  0.885  aQ  Z 


where  aQ  *0.5292  x  10  cm  is  the  Bohr  radius.  Thus 


and 


.  2  2 
4a  p 


4.44  x  10"6  Z^3 

T(T+1) 


.  /  /nUl  .  ,  -,2  2  r  1  +4T(T  +  1)1 

nQ  -  (n/8)N  fdZr  L - - - 5—] 


°  T2(T+1)2 


(20) 


r  ,  ^  8.88  x  lO'6  Z^3  Q  n-l  H.  ^  8.88  x  10^ Z^"  1-1 

M+ - T(TTT1 - COS0max  ]  -[2+  - T (T+T) - j 


-672/3. 


Since  the  use  of  small -angle  scattering  theories  requires  the  assumption  that  no 
electrons  are  scattered  at  angles  greater  than  9^*/  the  "theoretical "  error  can  be 
estimated  by  the  ratio  nc>/q>. 

As  an  example7  Equation  20  is  used  to  determine  the  errors  fora  Monte  Carlo 
calculation  of  1000  electron  histories  at  1-Mev  energy  with  an  aluminum  target.  The 
factor  8. 88x10“^  Z^0/T(T  +  1)  can  be  neglected  in  Equation  20  so  that 


/  /  /,*  2  KI  ,  ^2  r  1  +  4T(T+1)  , 

neA-(ir/8)r0  NdZ  t - 5 

r(T+ir 


to  - coi8max)-^  -  1/2 J 
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-9  22  3 

For  Z  =  13,  the  screening  length  is  1.99  x  10  cm  and  N  =  6. 03  x  10  atoms/cm  , 

The  value  of  O^x  is  taken  as  30  degrees,  for  which  the  Moliere  theory  is  still  ap¬ 
proximately  valid.  The  average  electron  energy  is  used  for  each  increment.  The 
results  with  the  accompanying  statistical  errors  are  shown  in  Table  4. 

Although  the  example  shown  in  the  table  has  been  carried  out  to  extremely  large 
error  values,  the  results  that  are  appreciably  larger  than  10  percent  cannot  be  con¬ 
sidered  as  anything  more  than  a  rough  Indication  of  the  dependence  of  the  error  on  the 
electron  energy  because  the  assumption  that  no  energy  loss  occurred  within  an  incre¬ 
ment  was  made  in  writing  Equation  19.  However,  in  the  example,  the  increment 
thickness  is  already  one-tenth  of  the  electron  range  after  the  tenth  slab  so  that  the 
effects  of  energy  loss  are  no  longer  negligible.  Also,  the  average  energy  used  in  com¬ 
puting  the  theoretical  error  contains  an  uncertainty  due  to  the  statistical  error  in  the 
number  of  transmitted  electrons. 

Nevertheless,  the  monotonic  increase  of  the  theoretical  error  indicates  that  it 
may  be  necessary  to  take  successively  thinner  increments  as  the  electron  energy 
decreases.  A  Monte  Carlo  calculation  can  decrease  the  Increment  by  taking  the 
thickness  In  proportion  to  the  energy,  or  by  specifying  that  each  increment  correspond 
to  a  certain  fractional  loss  In  the  electron  energy.  Either  procedure  will,  of  course, 
require  considerably  more  computation  than  the  use  of  a  constant  thickness.  For 
example,  Berger  (Reference  54),  who  used  a  logarithmic  energy-dependent  spacing, 
took  nearly  50  steps  to  reach  the  energy  value  corresponding  to  the  seventeenth  slab 
of  the  table.  The  corresponding  theoretical  error,  however,  would  have  been  only 
five  percent. 

5.2  RADIATION  TRANSPORT 

The  data  generated  from  the  Monte  Carlo  program  and  existing  theory  were  used 
to  formulate  mathematical  models  for  the  transport  of  electrons  and  bremsstrahlung 
through  multilayered  shields.  A  series  of  three  electron  models  are  presented: 

(1)  the  Monte  Carlo  program,  (2)  a  dose  calculation,  and  (3)  an  energy,  LET  spectra, 
and  dose  calculation.  There  Is  significant  difference  in  computing  time  so  that  each 
model  offers  certain  advantages  in  accuracy  and  detail  as  a  function  of  computing 
time.  The  bremsstrahl ung  model  includes  a  volume  distributed  source  and  the  analysis 
of  multilayered  shields. 

5.2.1  Primary  Electrons 

5.2. 1. 1  Transmission 

The  dimensionless  length  and  energy  parameters  of  Seliger  (Reference  60),  and 
Makhov  (Reference  61)  were  used  to  obtain  a  mathematical  model  for  the  number  of 
electrons  of  a  given  energy  that  penetrate  a  given  shield  thickness.  In  the  initial 
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Table  4 


THEORETICAL  AND  STATISTICAL  ERRORS  FOR  TARGET  INCREMENTS 
1-Mev  Electron*  Normally  Incident  on  Aluminum. 

Increment  Thlckneu  2*5  x  10"2  g/cm^  (9.25  x  10”3  cm). 


Slab 

Theoretical  Error  -  % 

Statistical  Error  -% 

1 

4.60 

3.16 

2 

4.65 

3.17 

3 

5.43 

3.20 

4 

5.73 

3.25 

5 

6.31 

3.34 

6 

7.52 

3.45 

7 

8.35 

3.63 

8 

10.15 

3.83 

9 

11.05 

4.12 

10 

13.90 

4.50 

11 

15.35 

5.13 

12 

18.05 

5.75 

13 

22.70 

6.65 

14 

25.95 

7.90 

15 

36.30 

10.18 

16 

42.40 

15.10 

17 

73.50 

25.80 
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electron-transmission  code  study  (Reference  62)  the  expression 


T(x,  E,  Z)  =  exp 


(•[ 


E(0. 585  Z 


Vx) 


-0. 271  /  "7(z  "  25) 


,-0.24 


was  derived  based  on  the  dimensionless  model  for  the  fraction  of  electrons  of  initial 
energy  E  which  penetrated  a  shield  of  material  Z  and  thickness  x.  This  expression 
was  derived  from  Monte  Carlo  data  for  beryllium,  aluminum,  iron,  and  silver  using 
1-,  3-,  4-,  6-,  and  8-Mev  normally  Incident  electrons.  The  dose  from  a  spectrum  of 
Incident  electron  energies  $e(E)  is  expressed  as 


De  -  F  j  K<|>e(E)T(E)^(E)dE 

where  K  converts  Mev/gm  to  rad  and  dE/dX(E)  is  the  energy  loss  in  the  target  for  an 
electron  of  energy  E.  The  limitations  of  this  expression  are: 


1)  The  shield  is  restricted  to  one  material; 


2)  The  electrons  are  assumed  to  penetrate  the  shield  and  retain  their  initial 
energy  for  dose  deposition;  and 

3)  The  electrons  are  assumed  to  be  incident  normal  ly,  but  are  corrected  for  num¬ 
ber  crossing  the  vehicle  surface  (F  factor). 


The  first  phase  of  the  electron  penetration  section  of  this  study  was  to  develop  a 
method  to  estimate  the  number  transmission  (ratio  of  penetrating  electrons  to  incident 
electrons)  through  multilayered  shields.  Approximately  2  hours  of  computer  time  was 
used  to  obtain  Monte  Carlo  penetration  data  for  monoenergetic  and  fission-spectrum 
electrons  having  normal  incidence  on  various  combinations  of  shield  materials. 

Analysis  of  these  data  indicated  that  the  number  transmission  for  the  total  shield  can 
be  obtained  by  the  following  scheme: 

1)  The  number  penetrating  the  first  slab  is  computed  with  the  "single  slab" 
method; 

2)  The  penetration  through  the  material  of  the  second  slab  is  also  computed  by  the 
single  slab  method  using  an  "effective  thickness*  " 

3)  The  penetration  through  any  slab  I  of  material  j  is  computed  using  an  effec¬ 
tive  thickness  equal  to  the  thickness  of  slab  t  plus  an  equivalent  thickness  for  slab  t, 
which  gives  an  attentuation  equal  to  that  of  slabs  1  through  1-1. 

4)  The  number  transmission  of  electrons  penetrating  a  multislab  shield  composed 
of  n  layers  of  materials  of  total  thickness 
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n 

x„  -  IE  *i 

1-1 

If  given  by 

T(x,,  Z,  |  1-1,2— n;  Ec)  -  T^,  Zn) 

where  by  definition  X  -  xn  +  .  The  Indlvlduol  thickness  of  the  n^  slab  Is  xn 

and  5^  It  an  equivalent  thickness  for  material  n  which  will  yield  a  transmission 
equal  to  that  penetrating  the  first  n-1  slabs  of  the  multi-layered  shield.  The  quan- 
tlty  Is  determined  from  the  Implicit  equation 

T<V,,  Eo'  ^  "  T(x,|  1-1,2,  — «-1,  Zn.„  E0) 

where  the  right-hand  side  of  the  equation  is  determined  from  the  previous  iteration. 
Equation  21  is  used  to  evaluate  the  left-hand  side  iteratively  until  a  value  of  is 
found  which  satisfies  the  equation. 

In  each  case  the  thickness  Is  in  terms  of  an  areal  density  in  units  of  grri/cm^.  This 
method  assumes  that  an  equivalent  thickness  of  any  material  that  produces  an  equal 
attentuatlon  of  electrons  will  also  produce  an  identical  energy  and  angle  distribution 
of  electrons  at  that  point.  Obviously  this  Is  not  exact  since  the  ratio  of  radiation  to 
Ionization  loss  is  dependent  upon  atomic  number.  However,  for  an  estimate  of  the 
number  transmission,  no  better  approximation  could  be  developed.  Results  of  the 
second  phase  of  the  study  revealed  an  alternate  method  of  comparable  accuracy.  A 
detailed  discussion  of  that  method  will  be  presented  In  the  next  section. 

The  accuracy  of  the  slngle-elab  expression  was  found  to  be  within  10  percent  In 
the  range  of  Z(  from  l  -  10  to  50.  Table  5  compares  the  results  of  the  multislab 
analysis  using  the  empirical  method  with  Monte  Carlo  results.  Thereare  no  experimen¬ 
tal  data  available  for  comparison  with  these  data*  Figure  34  Illustrates  the  improve¬ 
ment  obtained  using  the  multislab  method  compared  to  the  single-slab  method.  The 
transmission  for  a  multislab  shield  will  be  bracketed  by  the  single -slab  transmission 
curves.  For  this  reason,  there  is  no  advantage  to  using  a  multislab  approximation  when 
the  atomic  numbers  of  the  shields  are  similar  since  both  single  transmission  curves  are 
close  to  the  multislab  case.  A  shield  of  beryllium  and  lead  would  be  maximum  tests  of 
cases  where  the  single  slab  Is  not  representative  of  the  multislab  transmission.  Using 
the  lighter  atomic  weight  material,  the  maximum  error  in  the  single-«lab  method  Is  a 
factor  of  two,  while  use  of  the  heavier  atomic  weight  material  produces  factors  of  12. 
The  multislab  approximation  is  in  error  by  less  than  a  factor  of  two,  but  at  large  frac¬ 
tions  of  the  total  range  the  errors  again  become  large. 

It  Is  concluded  that  for  shield  thicknesses  less  than  one-half  the  electron  range, 
the  multislab  approximation  provides  15  percent  accuracy  In  the  calculated  number 
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Table  5 

MULTI  SLAB  ELECTRON  NUMBER  TRANSMISSION 


E 

Slab 

No. 

Material 

Thickness 

Total  Thickness 

Transmission 
(Monte  Carlo) 

f 

1 

Be 

.51  gm/cm^ 

.51  grr/cm^ 

.358 

.378 

1 

2 

Pb 

.51 

1.02 

.092 

.052 

s 

3 

Be 

.51 

1.53 

.020 

.026 

s 

4 

Pb 

.51 

2.04 

.003 

.003 

I 

5 

Be 

.51 

2.55 

.0002 

.001 

o 

n 

f 

1 

Pb 

.51 

.51 

.107 

.134 

I 

2 

Be 

.51 

1.02 

.03 

.065 

s 

3 

Pb 

.51 

1.53 

.01 

.019 

s 

4 

Be 

.51 

2.04 

.002 

.008 

I 

5 

Pb 

.51 

2.55 

— 

■  ■■■ 

o 

n 

f 

1 

Al 

.17 

.17 

.517 

.531 

l 

2 

A! 

.17 

.34 

.363 

.326 

s 

3 

Al 

.17 

.51 

.236 

.221 

s 

4 

Al 

.17 

.68 

.156 

.153 

I 

5 

Al 

.17 

.85 

.103 

.099 

0 

n 

*  Using  multislab  method. 


Transmission 


FIGURE  34.  Comparison  of  single  and  multislab  results  for 
Be-Pb-Be-Pb-Be 
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transmission.  The  maximum  errors  occur  at  large  percentages  of  the  electron  range 
when  there  are  large  differences  in  the  atomic  numbers  of  the  shield  materials.  There 
is  no  simple  modification  that  can  be  made  to  improve  this  method.  However/  an 
alternative  method  is  discussed  later. 

The  second  phase  of  this  study  was  the  development  of  a  mathematical  expression 
for  the  energy  spectrum  of  the  penetrating  electrons.  Monte  Carlo  electron- 
penetration  runs  were  made  for  beryllium/  aluminum,  iron,  silver,  and  lead  for  elec¬ 
tron  energies  from  0.4  to  6  Mev.  Both  normal  incidence  and  isotropic  distributions 
were  studied. 

The  derivation  of  an  empirical  expression  to  describe  the  differential  number 
spectrum  of  electrons  penetrating  a  shield  required  a  large  effort.  A  chronological 
summary  of  the  various  methods  will  be  presented  to  demonstrate  the  possible  mathe¬ 
matical  models  and  their  restrictions. 

The  first  method  was  based  on  the  observation  that  the  differential  number  spectra 
for  a  given  energy  group  appeared  to  be  a  family  of  Gaussian  curves.  As  the  final 
energy  decreased,  the  maximum  of  the  Gaussian  decreased.  Its  center  shifted  to  a 
larger  thickness,  and  its  half-width  at  half-maximum  Increased.  A  typical  family  of 
curves  for  a  1-Mev  electron  normally  incident  on  aluminum  is  shown  in  Figure  35. 

For  aluminum  shields  this  family  of  curves  was  represented  by  the  expression 

3 

♦,(E,x)  "  .xp C -2. 2  x  103.2-35  (x  -4.5  +  0.41.)2] 

where  c  is  dimensionless  energy,  the  ratio  of.  the  final  to  initial  energy,  and  x  I* 
dimensionless  length,  the  ratio  of  the  thickness  x  to  the  range  of  the  initial -energy 
electron.  As  this  analysis  was  extended  to  other  materials,  the  Gaussian  curves  began 
to  skew.  Attempts  to  correlate  data  using  a  translation  plus  a  rotation  of  the  skewed 
Gaussian  yielded  an  expression  too  complicated  for  practical  application.  This  ap¬ 
proach  was  subsequently  abandoned,  even  though  a  good  representation  was  achieved 
for  aluminum  shields. 

The  second  method  was  based  on  the  shape  of  the  differential  number  spectrum 
plotted  for  various  aluminum  shield  thicknesses.  The  Monte  Carlo  data  provided  a 
histogram  such  as  shown  In  Figure  36.  Experimental  data  of  White  and  Millington 
(Reference  63)  were  similar  to  these  data  and  are  shown  in  Ftgure  37. 

Since  the  quantities  of  Interest  to  be  computed  are  dose  and  LET,  several  assump¬ 
tions  were  made  In  obtaining  a  mathematical  model  for  the  differential  number  and  LET 
spectra.  Hie  first  assumption  was  that  a  curve  of  the  form  shown  in  the  previous  two 
Figures  (36  and  37)  can  be  represented  by  an  exponential  as  shown  In  Figure  38.  The 
approximation  Is  defined  as  zero  in  the  region  Emax  to  E0,  where  Emax  is  defined  as 
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aluminum  for  1-Mev  Incident  electrons 


0  0.2  0.4  0.6  0.8  1.0 

Energy  (Mev) 


FIGURE  36.  Mont*  Carlo  histogram  for  1-M*v  electrons 
incident  on  aluminum 
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FIGURE  37.  Transmission  through  mica  of  the  monoenetgetic 
electrons  of  0. 2065  Mev  (Reference  63) 
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N(E)  (arbitrary  units) 


FIGURE  38.  Mathematical  model  for  differential  number  spectrum 


the  maximum  energy  an  electron  of  Initial  energy  EQ  can  retain  when  penetrating 
straight  forward  through  the  shield  thickness  .  By  adjusting  the  integral  under  this 
curve  so  that  the  integral  of  the  fit  is  equal  to  the  integral  of  the  Monte  Carlo  data, 
the  proper  number  LET  is  retained.  It  can  be  observed  from  Figure  39  that  the  LET  of 
the  higher-energy  electrons  is  a  constant,  and  the  assignment  of  an  electron  to  a 
slightly  lower-energy  bin  does  not  effect  its  LET.  The  other  constraint  is  that  the 
integral  of  the  total  energy  in  both  curves  should  be  equal  to  retain  the  same  energy  or 
dose.  Both  restrictions  cannot  be  satisfied  simultaneously;  however,  a  compromise 
that  created  less  than  5  percent  error  was  obtained  for  aluminum  shields. 

The  procedure  for  this  analysis  was  to  analyze  the  histogram  for  the  number  of 
electrons  in  each  energy  bin  penetrating  a  shield  as  a  function  of  initial  energy,  shield 
material,  and  thickness.  A  value  for  the  differentia!  number  spectra  at  the  mid  point 
of  each  bin  was  obtained  by  dividing  the  number  of  electrons  in  a  bin  by  the  bin  width. 
The  differential  number  spectra  was  assumed  to  be  represented  by  the  relation 

♦.(E)  -  A.xp[-b(E0-E)]  E  <  Emax 

♦  (E)  -  O  E  >  Em.,, 

o 

where  the  constants  A  and  b  are  dependent  on  the  atomic  number  of  the  shield,  the 
initial  electron  energy,  and  the  thickness  of  the  shield.  Values  for  b  were  obtained 
from  the  shape  of  $e(E)  plotted  versus  EQ-E.  The  constants  A  were  obtained  by 
averaging  the  integral  of  the  expression  over  the  limits  of  each  bin  to  yield  the  number 
and  energy  in  each  bin.  The  Initial  analysis  provided  the  following  relations  for  A 
and  b: 

A  =  1.6(1  -0. 108  E0)  (x/E0)2*5  exp(  -0. 0435  Z) 

(VE0)_1* 25  0. 25  -  8.5  x  10"3 Z)/E0  x  >  0.325  E0 

b  *  < 

4.1(1.25 -8.5  x  10'3Z)/Eo  x_*  0.325  E0 

This  set  of  constants  was  satisfactory  for  small  x/tQ  but  yielded  laige  values  at  higher 
x/EQ.  Also  the  values  were  too  small  at  high  Z's. 

Re  -evaluating  the  data,  the  large  values  at  high  x/iQ  were  corrected  by  employ¬ 
ing  a  double  exponential  fit  for  A  which  corrected  the  error  at  large  */EQ.  The 
revised  constants  were 
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FIGURE  39.  LET  spectrum  of  various  radiation  types  (Reference  67) 


A  =  4.4  x  TO3  z"^*2^  Eq_1  * 23  exp(  -141  x/EQ)  +  1. 35  x  105  exp(-O.06  Z) 
xEo"°* 74  ®XP(-21  X/E0) 

b  -  (x/E,,)-1*46  (1.53  -0.0104Z)/Eo  (22) 

This  set  of  constants  was  cumbersome  and  the  expression  for  differential  spectra  was  too 
sensitive  to  changes  in  The  b  values  for  small  x/EQ  were  large  numbers,  and 

a  small  error  in  E^x  could  affect  the  exponent  markedly.  The  final  expressions  em¬ 
ployed  subroutines  to  compute  A  where  the  total  number  penetrating  the  shield  was 
determined  by  the  original  number  transmission  relation  Equation  21  and  then  A  com¬ 
puted  from  the  relation 

_ Ntrans  b _ 

A  exp(-bE0)  [exp(bEmax)  -  1] 

In  each  case  b  was  found  to  be  represented  by  Equation  22. 

The  results  were  found  to  be  in  agreement  with  the  Monte  Carlo  calculations. 

5. 2. 1.2  Results 

Calculations  using  fission  spectrum  electrons  and  aluminum  shields  produced  results 
within  15  percent  agreement  of  the  Monte  Carlo  calculations.  This  method  was  ex¬ 
tended  to  muitislab  shields  with  better  agreement  than  the  previous  methods.  There 
remains  one  assumption  in  the  mathematical  model  which  tends  to  overestimate  the 
final  dose.  All  electron  transmission  data  were  fitted  to  normal -incidence  penetration 
data.  Figure  40  shows  the  number  transmission  of  l-Mev  electrons  through  aluminum 
for  0,  30,  and  60  degrees  and  isotropic  incidence.  (0  degree  is  normal  incidence.) 

The  Isotropic  number  transmission  was  found  to  be  about  50  percent  of  the  normal 
transmission,  but  the  equilibrium  spectrum  was  achieved  with  a  smaller  thickness.  This 
softer  spectrum  caused  the  dose  per  transmitted  electron  to  be  greater  for  isotropic  dis¬ 
tribution  than  for  normal  incidence.  Comparing  the  dose  transmission  from  normal  and 
isotropic  incident  fluxes  the  values  were  similar  and  the  larger  normal  transmission 
factors  were  used  to  retain  a  safety  factor.  The  integration  of  the  normal -transmission 
factors  over  the  slant  paths  encountered  by  an  isotropic  distribution  was  processed  for 
l-Mev  electrons  on  aluminum.  This  method  produced  the  same  transmission  for  iso¬ 
tropic  flux  as  the  Monte  Carlo;  however,  the  multiple-scattering  theory  is  so  uncertain 
and  the  calculations  are  so  lengthy  that  this  method  was  not  programmed. 
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5. 2. 1.3  Dose  Calculation 

The  electron  dose  transmission  factor  De(n)  for  a  given  sector  composed  of  n  layers 
is  expressed  as 

D.w  =  K  $  ST®  j*.(E„)Tn<Eo-  E)dEodE 

where  K  =  conversion  factor  Mev/gm  to  rad 

dE  -1  2 

-r-(E)  =  rate  of  energy  deposition  for  electron  of  energy  E,  Mev  gm  cm 

ax 

<|>e^0)  =  spectrum  of  incident  electrons  normalized  to  one 

T  (E0-*.  E)  =  transmission  of  electrons  with  energy  E  through  shield  (obtained  by 
use  of  Equation  21,  Section  5.  2. 1. 1) 

This  operation  (1)  computes  the  number  of  electrons  of  energy  E  penetrating  the  shield 
from  incident  electrons  of  energy  Eg,  (2)  integrates  over  all  incident  energies  to  ob¬ 
tain  the  total  number  of  electrons  of  energy  E  penetrating  the  sector,  (3)  computes 
the  energy  deposited  by  these  electrons,  and  (4)  integrates  over  all  final  electron 
energies  to  compute  the  total  dose. 

The  total  dose  transmission  through  all  sectors  is  expressed  as 

D#  =  XZ  D0(n)  S#  F^(n)  -A_  (n) 
n 

where  D#(n)  is  defined  above,  Se  Is  the  free -space  flux,  Fe(n)  is  the  fraction  of  the 
free-space  flux  incident  on  sector  n,  _A_(n)  is  the  fraction  of  electrons  penetrating 
sector  n  which  are  incident  on  the  dose  point,  and  up  is  the  fractional  solid  angle. 

The  function  F#(n)  has  been  shown  by  Evans  (Reference  64)  to  be  one-fourth  for  an 
isotropic  distribution.  An  analysis  of  the  value  of  F0(n)  for  other  distributions  will  be 
presented  in  the  section  7.0  discussing  the  effects  of  the  incident  angular  distribution. 

It  has  previously  been  shown  that  for  electrons  a  diffusion  equilibrium  is  obtained  where 
the  angular  distribution  can  be  described  by  cos^O  (9  is  the  scattering  angle  from  the 
normal  to  the  shield).  The  function  J\_(n)  was  previously  discussed  as  a  function  of 
cos^9  in  Section  5. 1.3  as 

_A_  (n)  ■  fQ  cos^Q 
where  fQ  is  a  normalizing  factor. 
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The  Monte  Carlo  program  for  electron  penetration  ha*  been  modified  for  use  in 
computing  a  detailed  analysis  of  electron  penetration  in  a  vehicle.  All  the  output 
such  as  number  spectrum/  LET  spectrum,  energy  loss  per  unit  path  length,  and  dose  can 
be  computed  with  the  Monte  Carlo  electron  code  as  well  as  the  simpler  code  using  the 
empirical  fit.  For  single  material  or  shields  where  there  is  only  one  sector,  the  Monte 
Carlo  will  provide  a  much  more  accurate  answer  in  the  same  computing  time. 

5. 2. 1.4  LET  Calculation 

The  calculation  of  the  energy  imparted  by  each  primary  electron  per  unit  path 
length  in  a  given  medium  to  obtain  the  Linear  Energy  Transfer  (LET),  was  based  on  the 
mean  energy  loss  values  of  Nelms  (Reference  51).  An  estimate  of  the  secondary,  ter¬ 
tiary,  and  higher-order  products  was  made,  but  there  are  insufficient  data  to  justify 
accurate  calculations.  It  is  believed  that  the  small  correction  which  might  be  obtained 
would  not  justify  the  time  required  for  the  calculation.  In  calculations  of  permissable 
exposure  levels  the  1963  report  of  the  Relative  Biological  Effect  Committee  (Reference 
65)  recommends  the  use  of  the  concept  of  the  Quality  Factor  QF  in  terms  of  a  LET., 
defined  for  primary  particles.  It  indicates  that  for  X,  y,  and  13  rays,  the  QF  is  close 
to  unity  and  the  LET  is  always  low  except  for  very  low  electron  energies.  The  concept 
that  X,  y,  and  13  radiation  has  a  low  QF  or  RBE  has  been  recognized.  Cormack  and 
Johns  (Reference  66)  and  Burch  (Reference  67)  have  made  detailed  calculations  of  LET 
spectra  to  verify  this  point.  It  was  noted  by  Schaeffer  (Reference  68),  who  compared 
the  proton  LET  spectrum  with  the  X-ray  LIT  spectrum,  that  the  spectrum  for  LET  greater 
than  10  kev/micron  is  an  artifact  for  X-rays  since  the  mechanism  for  energy  dissipation 
below  1  kev  is  incompletely  understood.  The  maximum  LET  for  electrons  has  been  esti¬ 
mated  to  be  38  kev/micron,  and  a  limiting  value  of  30  kev/micron  has  been  recommen¬ 
ded  by  Schaeffer  as  the  dividing  line  between  high  and  low  LET. 

Figure  39  is  a  representation  of  the  data  of  Burch  (Reference  67)  smoothed  over 
energy.  Since  electrons  deposit  their  energy  in  much  the  same  manner  as  X-rays  (i.e. , 
by  production  of  low-energy  electrons)  the  QF  of  electrons  (Reference  68)  should  be 
that  of  X-rays,  namely  unity.  Figure  41  shows  a  comparison  of  the  LET  spectra  of  X- 
rays,  solar  protons,  and  recoil  protons. 

5.2.2  Bremsstrahlung  Shielding 

5. 2. 2.1  Transmission 

The  original  computer  program  developed  at  Boeing  to  compute  the  bremsstrahlung 
dose  in  a  space  vehicle  (Reference  62)  had  two  assumptions  which  limited  the  accuracy 
of  the  calculation  to  an  order  of  magnitude  in  some  situations.  These  assumptions  were 
(1)  the  photons  were  ail  produced  on  the  incident  surface  of  the  shield,  and  (2)  an  in¬ 
finite  slab  thickness  was  used  to  predict  the  photon  yield.  The  first  assumption  is  in 
error  when  a  shield  Is  composed  of  a  thin  low-atomic-number  shield  followed  by  a 
high -atomic -number  shield.  Since  the  bremsstrahlung  production  Is  proportional  to 
atomic  number,  a  low  production  is  predicted  even  though  electrons  can  penetrate  the 
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FIGURE  41.  LET  spectrum  of  X  rays,  solar  protons,  and  recoil  protons 
from  neutrons  (Reference  68) 
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first  layer  and  produce  more  photons  in  the  second  layer.  The  second  assumption  pro¬ 
duces  errors  in  thin  shields  since  only  part  of  the  electron  energy  Is  deposited  in  the 
shield.  Both  these  assumptions  were  eliminated  by  the  modifications  made  during  this 
study. 

The  original  codes  employed  the  following  relation  for  bremsstrahlung  dose 


[  4(1 


ln(j^-)3  dE  d(hV) 


where  Dg  is  the  bremsstrahlung  dose  created  in  an  N-layered  shield  by  an  incident 
spectrum  of  electrons  <|>_(E).  K(hv')  is  the  flux-to-dose  conversion  factor  for  photons 
of  energy  hY”,  Bj(pj(hY)  Xj)  is  the  buildup  factor  for  photons  of  energy  hV  for 
pxj  mean  free  paths  in  shield  i,  pj  Is  the  mass  attentuatlon  coefficient,  and  x|  is  the 
slab  thickness.  The  expression 


W(E) 


1.98  x  104(1. 96  E  +  2) 


ZE 


(24) 


is  the  average  amount  of  energy  radiated  by  bremsstrahlung  in  a  thick  target  as  shown 
by  C.  S.  Wu  (Reference  69).  The  term  in  the  brackets  is  the  energy  radiated  per  unit 
energy  interval  derived  by  Wyard  (Reference  70). 


I(hr,  E)  =  t  4  (1  -  -jj—  )  -  ^-3  ln(E/hV)  (25) 

The  electron  Monte  Carlo  code  was  used  to  obtain  data  for  bremsstrahlung  source 
distribution.  Data  for  aluminum,  iron,  silver,  and  lead  for  1-,  3-,  and  6 -Mev  electrons 
were  obtained.  It  was  generally  observed  that  the  rate  of  bremsstrahlung  production 
was  approximately  constant  over  the  region  where  the  bremsstrahlung  of  a  particular 
energy  group  could  be  produced. 
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The  volume-distributed  bremsstrahlung  was  computed  by  modifying  the  original 
expression  (Equation  23).  Consider  a  spectrum  of  incident  electrons  and  select  the 
number  of  electrons  of  energy  E  per  unit  energy  interval,  <f>e(E).  From  Equation  24 
tor  the  energy  radiation  per  electron  W(  E,  Z)  in  a  thick  target,  the  photon  energy 

release  is 


photon  energy  release 

from  =  $e(E)  W(E,  Z) 

electrons  of  energy  E 

The  fraction  of  these  releases  as  photons  of  energy  hV  is  given  by  Equation  25 
defined  as  l(hv,  E).  Thus,  the  photon  energy  of  hy  released  by  electrons  of  energy 
E  is 


photon  energy  release 

of  energy  hV  from  *  I  (hy,  E)  <|>e(E)  W(E,  Z) 

electrons  of  energy  E 


Since  only  electrons  with  energy  greater  than  hV  can  create  photons  of  energy 
hV~,  the  thickness  XQ(E,  hV)  where  photons  can  be  created  can  be  defined  by 


hy  -  E0  - 


A 

I 


3x  (E©'  x>  dx 


where  dE/cb^E^  x)  is  the  rate  of  energy  loss  of  an  electron  of  Incident  energy  E0  at  a 
thickness  x  in  the  shield.  For  this  calculation,  the  electron  path  is  normal  to  the 
shield  to  obtain  maximum  penetration. 


The  source  strength  per  unit  volume  of  photons  of  energy  hy  produced  by  elec¬ 
trons  of  energy  E  is  then 


l(hvyE)  »,(E,x)  W(E,  Z) 
X0(E,hV) 


Hx/Xq) 


where  F(x/XQ)  can  be  a  distribution  of  the  source  along  x. 

The  transmission  of  the  radiation  through  the  source  slab  n  and  through  the  re¬ 
maining  slabs  is  defined  as 

N 

T(hV,  Xyy  )  Xj ) 
i  ■  n+1 
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where  h  V  is  the  energy  of  the  penetrating  photon,  xn  i*  the  thickness  of  the  radiating 
slab  and 

N 

xi 

i  =nH 


is  the  thickness  of  the  remaining  slabs.  Thus,  the  photon  energy  reaching  the  receiver 
from  photons  of  energy  hV  created  by  electrons  of  energy  E  in  slab  n  is: 


<t>(E,hZ)  = 


N 

l(hV,E)<t>e(E)  W(E,Zn)  T(h\f,xn  YL  *i 
_ t  »  n+1 

X0n(E,  hV) 


Integrating  over  all  incident  electron  energies  will  yield  the  total  photon  energy 
reaching  the  receiver. 

t(hV)  -  ^(E,  hV)dE 

Multiplying  this  integral  by  the  dose  conversion  factor  for  energy  h/",  K(hV)  and  then 
integrating  over  all  photon  energies  gives  the  total  dose  from  all  photon  energies 
created  by  <|>e(E)  as 

Db  =  J  K(hV)  A(hv)  *{bsr)  d(h\f) 

where  A(hV)  is  a  distribution  function  for  angle  as  a  function  of  photon  energy. 

The  $e(E)  penetrating  slab  n  end  incident  on  slab  n  +  1  Is  computed  using  the 
electron  code  (Equation  21).  This  n+j(E)  •*  u,ec*  to  compute  the  dose  D^+j  and 

the  computation  continued  until  the  dose  from  photons  created  In  each  slab  is  computed. 

Each  sector  is  assumed  to  be  composed  of  increments  of  solid  angle  0.05  stera- 
dians.  The  dose  from  a  truncated -cone  volume  source  composed  of  N  layers  has  been 
derived  by  Foderaro  and  Obenshain  (Reference  71).  For  a  volume -distributed  source 
where  the  source  strength  is  constant,  the  function 

N 

T(hZ,Xn,  n  x.) is  given  by 

j=n+1  J 
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N 


Tfl*,Xn, 


j*n+l 


EA>_sec0) 

Zpn(hr)  C  E2(b^  "  E2(b3)  +  seed  ' 


SV  B(bJ 
n  o 


where  Pn(hv)  =  attenuation  coefficient  for  photon  of  energy  in  slab  n 
Xn  =  thickness  of  slab  n  containing  the  source 
0  =  subtended  angle  (see  Figure  42) 


N 


H.(hY-)  X. 


N 


p.(hV)  X. 


i=*n 


r<*> 

E2(b)  *  b  ^  exp(-t)/t2  dt 


B(b^)  *  buildup  factor  =  A^exp  -  A^  +  A2exp  -  a^ 


a2  *  1  -  °1 


a^,  a2  as  defined  in  Reference  62. 


If  the  volume  source  is  a  linear  sum  of  exponentials,  such  as  Sy  (x)  =  5 

the  transmluion  is  defined  as  1 

T(hv)  -  L.-ypp  SI 


where 


8j  -  exp(kjxn)  CEjfcj)  -  E^secS)] 


p  +k.  p  +k. 

-  exp^bj/pn)  t  E1  {\>]  — )  -  E1  (bg  -n^  -)  3 


2(b^sec8) 
seed  ^ 


j  exp(  kj  x) 
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-  E^bg)  -  E^(b3  sec9) 


+  exp(k.b3/Hn)  t^1(b1  — 


p  $ec9  +  k,  p  sec9  +  k. 

>-E,<b3-^ - 1)1 


E^b) 


exp(-t) 

t 

b 


Special  cases  when  (1)  bj  /*0,  kj  31  0;  (2)  bi  =0,  k.  f  0;  (3)  b.  /  0,  k.  =  0; 

(4)  b]  0,  kj  f  -p,;  (5)  b]  ?  0,  k.  =  -p,  sec9;  and'(6)  bl  f  0,  k|  f  0,  -p* 

~PjSec9  were  evaluated  and  are  summarized  in  Table  6.  These  special  cases  require 
additional  analysis  of  the  general  expression  since  otherwise  division  by  zero  would 
occur.  The  quantity  b]  is  the  number  of  mean  free  paths  of  shielding  following  the 
source  layer  and  kj  is  the  is  the  coefficient  of  the  equipment  of  the  i"1  term  in  the  series. 
Sy(x)  “ZlSj  exp(kjx).  The  quantity  B  is  the  argument  for  the  i*h  term  in  the  series. 

In  Case  1  (b]  /  0,  kj  =  0)  there  is  shielding  beyond  the  source  layer  but  the  i^  source 
term  is  unity.  Case  2  (bj  *  0,  kj  4  0)  has  no  shielding  after  the  source  layer  but  the 
i^1  source  term  has  a  finite  value.  Case  3  (b]  -  0,  kj  *  0)  has  no  shielding  and  the 
Ith  source  term  is  unity.  Case  4  (bj  0,  k;  =  -p#)  has  a  shield  but  the  coefficient  of 
’  the  i^1  source  term  is  equal  to  the  photonattenuation  coefficient  p,.  Case  5  (bj  *0, 
kj  *  -p|Sec9)  has  a  shield  but  the  coefficient  of  the  i™  source  term  is  equal  to  the 
photon  attenuation  coefficient  times  the  secant  of  the  subtended  angle.  Case  6 
(b^  *0,  kj  0,  f  -pj,  /  -p|Sec9)  has  a  shield  and  the  coefficient  of  the  i^  source 
term  does  not  have  the  values  0,  -p,,  or  -pssec9. 

The  angular  distribution  of  bremsstrahlung  has  been  described  in  the  mathema¬ 
tical  model  by  the  function  F(hv/  E,  9  =  0,  %)  which  represents  the  fraction  of  the 
total  bremsstrahlung  created  by  an  electron  of  energy  E  incident  at  an  angle  9Q  which 
will  have  an  energy  hv  and  be  directed  in  the  direction  of  9  *  0.  This  function  is 
required  for  a  thick  target,  since  the  differential  cross  sections  are  only  useful  in  a 
Monte  Carlo  or  transport -type  analysis.  The  existing  data  are  integrated  over  the 
final  energy  and  data  are  reported  as  energy  per  square  centimeter  (References  72,  73, 
and  74).  Typical  data  are  shown  in  Figure  43.  Analysis  of  these  data  indicates  that 
a  sum  of  exponentials  may  be  an  approximation  to  the  distribution  function.  There 
were  insufficient  data  to  develop  a  successful  expression;  however,  the  machine  pro¬ 
gram  has  the  flexibility  to  accept  any  analytical  expression  developed.  Until  more 
information  is  available,  a  cos*0  distribution  appears  to  be  a  reasonable  approximation 
to  the  bremsstrahlung  distribution.  If  the  angular  distribution  of  the  incident  radiation 
is  assumed  to  be  isotropic,  the  radiation  produced  will  also  be  Isotropic  over  2w  and 
independent  of  the  scattering  distribution.  This  situation  Is  usual,  but  an  analysis  of  a 
nonisotropic  incident  distribution  will  be  presented  In  a  separate  discussion. 
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Table  6 

SUMMARY  OF  SPECIAL  CASES 


Case  1 


Case  2 


Case  3 


Case  4 


Case  5 


b,  /  0,  k.  =  0 


Aj  E.(b3sec0  E,(b.sec0) 

>  - is 1  w  -  w +  -  J_U_ 


sec0 


b,  "  0,  k.  /  0 


A. 


r 


SiS^  j  exp(kfh)  In(sec0)  -  [  E^bg)  -  E^bg  sec0)  ] 


L 


b,  -  0,  kj  =0 


A.  E_(b.  sec0)  . 

1  “  2m,  C  1  ‘  E2(b3)  £35  £3&] 


2Ms 

bT  ^  k,  =  -ms 

A. 


Bi  "aCT  ®Xpkjh  C  El(bl)  "  El(bl  *ec0)]  "  t  Ef(b3)  -  E^bgsece)] 

+  •Xp(k?b3/M,)^  -  E1  C  bg(,ec0  -  1)]  +  E^b^cO  -  1)]  -  ln(£l)  "j 

bl  ^  0/  *  -H|Sec0 

Aj  f 

°!  '  'a^'}  •xPklhCEl<bl>  -  E,<b,  we6)l  -  [Ejttj)  -E^mcS)] 


klb3  f 
-  exp(-~-)  -  E. 

p* 


Cb3(l  -S6C0)]  +  Ej  Cb^l  -sec0)  ]  -  In  (bj/bg)  ^ 
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Table  6 
(Continued) 


Case  6  b 
G 


1  /  0,  kj  /  0,  /  li,  *®c0 


Ai 

1  2^ 


+  exi 


exp(k.h)  L  E^)  -  E](b1  sec0)3  -  exp(kjh)  [  E^bj)  -  E^secO)  3 


..  M  sec0+k. 

M‘‘ec9  -El[b3(A_ - — )3 


k.b,  C  Pssec0+k. 

p+k.  Ms+ktx-,l"l 

-  E,  t  hll-ir2'1  +  E1  1  b3(  M,  || 
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5. 2. 2. 2  Dose  Col culations 

The  mathematical  model  for  bremsstrahlung  shield  was  incorporated  in  the  origi¬ 
nal  vehicle  shielding  code  for  bremsstrahlung.  The  bremsstrah lung  dose,  Dg,  inside 
a  vehicle  is  computed  from  the  expression 

°B  *  TBi  A,  “I 

i 

where  $ej(E)  is  the  flux  of  electrons  incident  on  sector  i  of  the  vehicle.  Tg|  is  the 
bremsstrahlung  penetrating  sector  i  resulting  from  an  incident  flux  of  $e*(E)  electrons 
of  a  specific  spectrum  and  is  the  transmission  factor  discussed  previously.  Aj  is  the 
angular  weight  factor  which  is  the  ratio:  dose  scattered  to  the  receiver/isotropic  dose 
to  the  receiver,  and  u>  is  the  fractional  solid  angle  subtended  by  the  receiver.  A  dis¬ 
cussion  of  the  computation  of  A.  is  presented  in  Section  7,  which  describes  the  general 
effects  of  angular  distribution. 

5. 2. 2. 3  Results 

The  results  of  the  modified  bremsstrahlung  model  were  compared  with  the  data 
from  the  original  model.  These  data  are  shown  in  Table  7.  The  agreement  is  within 
five  percent.  Since  the  source  has  been  distributed  uniformly  for  each  energy  group 
of  photons,  the  amount  of  bremsstrahlung  penetrating  a  shield  should  be  slightly 
greater  than  the  surface  source.  The  difference  is  the  decrease  in  absorption  o^the 
photons.  For  an  aluminum  shield,  the  source  volume  is  approximately  3  gm/ cm 
thick.  This  would  reduce  the  1-Mev  photons  by  15  percent,  the  3-Mev  photons  by 
10  percent,  and  the  6-Mev  photons  by  8  percent.  For  shields  less  than  the  source 
volume,  only  a  fraction  of  the  total  bremsstrahlung  is  produced,  so  the  volume  source 
results  should  be  lower  than  the  surface  source.  The  exponential  series  distribution 
has  been  validated  with  hypothetical  constants;  however  no  meaningful  data  were 
generated.  Tables  8  to  11  show  results  obtained  from  this  code  for  lead,  iron,  alu¬ 
minum,  and  beryllium  shields  and  a  75-percentile  man  as  a  function  of  shell  thickness 
and  location  of  dose  point  in  the  man. 
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Table  7 

COMPARISON  OF  BREMSSTRAHLUNG  RESULTS 


Target 

Fission  Spectrum 

Fission  Spectrum 

All  Incident  Electrons  Stopped  on  Surface 

Uniform  Distribution* 

2 

1  gm/cm  Be 

3.44  x  10  ^  rads/e/cm^ 

-9  2 

2.3x10  rads/e/cm 

1  gm/cm^  Al 

1. 15  x  10-8 

°r 

o 

X 

-«r 

e 

O 

2 

5  gm/cm  Pb 

4. 75  x  10"8 

4.9  x  10"8 

*  Photons  of  energy  E  produced  over  a  depth-equivalent  electron  energy  loss 
of  E0  -  E  where  E0  is  Initial  electron  energy. 
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6.0  PROTON  SHIELDING 


Once  the  proton  flux  incident  on  a  space  vehicle  has  been  selected,  the  detailed 
interactions  of  ihe  incident  proton  flux  with  the  materials  of  the  space  vehicle  and  the 
astronaut's  body  can  be  considered.  The  objectives  of  this  phase  of  the  contract  were 
to  develop  techniques  for  calculating  the  absorbed  dose,  number-LET,  and  energy-LET 
spectra  at  arbitrary  points  in  a  space  vehicle.  A  method  for  calculating  semiconduc¬ 
tor  dislocation  densities  was  also  desired.  Two  computer  codes  were  developed  during 
this  study.  The  first  treats  those  situations  in  which  the  incident  proton  flux  is  the 
dominant  energy  depositor..  Basically,  the  existing  Boeing  computer  code  was  found 
to  be  suitable  for  this  application,  and  it  was  revised  to  generate  the  required  data. 

A  more  comprehensive  treatment  of  proton  interactions  was  included  in  a  second  com¬ 
puter  code,  which  was  developed  to  verify  both  the  approximations  made  in  the  first 
code  and  to  treat  situations  in  which  the  secondary  radiations  induced  by  the  incident 
proton  flux  are  significant. 

6. 1  INTERACTIONS  OF  PROTONS  WITH  MATTER 

As  protons  pass  through  an  absorber,  they  lose  energy  primarily  through  collisions 
with  the  electrons  and  nuclei  of  the  absorber's  atoms.  Proton  bremsstrahlung  is  a 
relatively  insignificant  energy-loss  mechanism  in  comparison  with  the  collision  pro¬ 
cesses  (Reference  75).  For  protons  of  low  energy,  say  10  Mev,  collisions  with  elec¬ 
trons  provide  the  dominant  energy-loss  mechanism.  As  protons  of  increasing  energy 
are  considered,  collision  processes  with  the  nuclei  of  the  absorber  assume  increasing 
importance.  This  includes  only  the  "first  order"  interactions.  Ultimately,  the  entire 
kinetic  energy  of  the  proton  will  be  dissipated  in  the  ionization  and  excitation  of 
electrons  by  charged  particles,  or  carried  away  from  the  point  of  interest  by  neutral 
particles.  It  will  be  shown  that  the  eneigy  spectrum  of  solar  flare  and  trapped  protons 
makes  collision  processes  with  electrons  the  primary  energy-loss  mechanism,  although 
for  thicker  shields  and  deep  body  points,  the  dose  from  secondary  radiations  is  not 
negligible. 

6. 1. 1  Interaction  With  Electrons 


Charged  particles  in  general— and  protons  in  particular— lose  energy  to  atomic 
electrons  in  collisions  that  are  conveniently  divided  into  two  classes,  hard  and  soft. 

A  quantity  H  (the  amount  of  energy  transferred)  defines  the  division  between  these  two 
classes.  Bethe  (Reference  76)  has  shown  that  the  energy  loss  per  unit  path  length 
for  "soft"  collisions,  Q  <  H,  is  given  by 


2n  N(Zel 
“FT" 

m  0  c 
o 


r 

2  2 

2m  Be  H 

In 

o 

(1  -  82)  I2 

—  — 

-  B 
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where 

Q  *=  energy  transfer  per  collision 

3 

N  *  number  of  electrons  per  cm  of  absorber 

m  =  rest  mass  of  the  electron 
o 

Ze  *=  charge  of  the  penetrating  particle 
I  =  the  mean  excitation  potential  of  the  medium 
H  =  the  maximum  energy  transfer  in  a  single  soft  collision 
G  *=  v/c 

v  *  velocity  of  the  incident  particle 
c  *  velocity  of  light 


For  H  sufficiently  large  that  the  binding  energy  of  the  electron  is  small  in  comparison, 
the  cross  section  for  an  energy  transfer  greater  than  H  is  essentially  that  for  a  free 
electron.  For  incident  proton  energies  less  than  93.8  Bev,  one  may  use  the  cross  sec¬ 
tion  for  spin  zero  particles,  and,  in  addition,  neglect  bremsstrahlung  energy  losses 
(Reference  77).  This  gives  an  energy  loss  per  unit  path  length  of 


o 


where 


Q 

max 


2  m  v 

0^7) 


from  relativistic  mechanics,  and  the  other  symbols  are  defined  as  before. 

Combining  the  two  expressions  gives  the  well-known  relation  for  the  total 
average  energy  loss  of  a  heavy  partich 
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It  is  convenient  in  shielding  calculations  to  use  density  independent  units  rather  than 
linear  thickness,  which  gives 


4ir  N  Z^e^ 
o 

T~ 2 

Am  B  e 


L 

r  2  2~ 
2  m  cT 
o 

- 

0  -  B2)  1 

J 

where 


23 

N  =  Avogadro's  number®  6.  025  (10  )  atoms/mole 
o 

A  «=  atomic  weight  of  the  absorber 
p  “  density  of  material 

In  what  follows,  thicknesses  will  be  given  in  the  density  independent  form  x. 

This  expression  has  been  found  accurate  if  the  following  conditions  are  ful¬ 
filled: 


1)  The  velocity  v  of  the  incident  particle  must  be  greater  than  the  largest 
electronic  velocity  of  the  absorber,  v  *  7.*  /\\  <  v.  For  oxygen,  this  means  that  the 
velocity  of  a  proton  must  be  greater  tlrcm  1  Mev.  Fortunately,  this  otherwise  quite 
serious  restriction  can  be  circumvented  by  the  use  of  a  correction  factor  c^  which 
takes  into  account  the  reduced  effectiveness  of  the  more  tightly  bound  electrons 
(Reference  78).  It  will  also  be  shown  that  in  space  vehicle  shielding,  the  importance 
of  the  low-energy  protons  is  small  if  a  few  gm/cm^  of  shielding  are  present. 

2)  The  velocity  of  the  particle  must  be  great  enough  that  its  effective  charge 
is  not  altered  by  the  capture  of  atomic  electrons.  There  ere  no  valid  theoretical 
expressions  for  (-dE/dx)  in  this  region,  and  experimental  results  must  be  used.  For¬ 
tunately,  this  energy  region  is  below  1  Mev  for  protons,  although  for  alpha  particles 
it  presents  a  more  serious  problem. 

Once  the  average  energy  loss  (-dE/dx)  has  been  acquired,  it  is  possible  to 
derive  range-energy  tables  in  the  manner  of  Rich  and  Madey  (Reference  79),  or 
Stemheimer  (Reference  80).  In  addition,  further  investigation  of  the  relative  impor¬ 
tance  of  the  hard  and  soft  collisions  gives  an  insight  into  the  exact  nature  of  energy 
dissipation  in  tissue.  By  taking  the  ratio  of  the  total  energy  loss  to  the  energy  lost 
in  soft  collisions,  we  find 


117 


By  examining  the  amount  of  energy  lost  by  collisions  of  greater  than  say  200  ev,  an 
insight  into  the  relative  amount  of  energy  dissipated  by  local  or  soft  collisions  can  be 
gained.  The  electrons  of  energy  greater  than  200  ev  have  been  shown  to  be  less 
effective  in  causing  biological  damage  than  the  soft  energy  loss  component  of  the  pro¬ 
tons.  Schaffer  has  examined  this  subject  in  detail  and  concludes  that  the  energy 
deposited  in  primary  collisions  should  be  divided  into  these  categories  for  an  accurate 
evaluation  of  radiation  effects  on  man  (Reference  81). 

6. 1.2  Interactions  with  the  Nucleus 

6.  1.2.1  Elastic  Collisions 

Protons  also  lose  energy  to  an  absorber  by  interactions  with  the  nuclei  of  the 
absorber's  atoms.  These  interactions  can  be  divided  into  elastic  and  inelastic  collisions. 
An  elastic  collision  is  defined  as  an  interaction  between  two  particles  in  which  the 
kinetic  energy  of  both  is  conserved  in  the  center-of-mass  coordinate  system.  It  can 
be  shown  that  the  energy  lost  by  elastic  scattering  due  to  Coulomb  or  nuclear  forces 
is  small  in  comparison  with  the  electron  ionization  energy  losses  (Reference  82). 

Only  in  the  case  of  the  light  nuclei  can  significant  energy  be  transferred  to  the  nucleus 
in  the  event  of  an  elastic  collision.  (The  heavier  nuclei  are  essentially  undisturbed 
by  the  passage  of  a  proton. ) 

One  basic  reason  which  renders  the  elastic  collision  processes  ineffectual  in 
transferring  energy  is  the  infrequency  of  large  angle  scattering  events.  At  high 
incident-proton  energies,  inelastic  nuclear  collisions  are  dominant  because  the  proton 
interacts  with  the  individual  nucleons.  At  low  energies,  where  elastic  collisions  are 
possible,  the  total  range  of  the  protons  is  so  small  that  only  a  small  fraction  of  the 
protons  have  elastic  collisions.  Coulomb  scattering,  which  is  possible  at  all  incident- 
proton  energies,  has  such  a  small  cross  section  for  large  angle  scatterings  that  it  is 
also  unimportant  as  an  energy  dissipation  process.  Only  in  the  case  of  hydrogen, 
where  elastic  collisions  alone  are  possible  up  to  the  pion  production  threshold  near 
300  Mev,  do  elastic  collisions  deserve  attention. 

The  low  probability  of  large-angle  elastic  scattering  events  has  another 
significant  result.  The  incident  protons  maintain  their  original  direction  while 
traversing  the  shield  to  a  very  close  approximation.  The  atomic  electrons,  with  their 
low  relative  mass,  do  not  deflect  the  incident  protons  through  large  angles  either  in 
single  collisions  or  by  multiple  scattering.  Multiple  scattering  resulting  from  Coulomb 
forces  is  also  small.  For  example,  the  mean  square  radial  spread  of  a  100-Mev  proton 
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beam  passing  through  2  gm/cm^  of  aluminum  is  0.01  cm  (Reference  83).  The  charac¬ 
teristic  of  protons  can  be  used  to  advantage  in  the  analysis  of  practical  shielding 
problems. 

6, 1.2.2  Inelastic  Collisions 

As  protons  of  higher  energy  are  considered,  the  inelastic  collision  process 
represents  an  increasingly  important  mechanism  for  modifying  the  incident  proton  flux. 
This  is  a  result  of  two  factors,  the  rapid  increase  of  proton  range  with  energy 
(R  *  KE^*^)  and  the  nearly  constant  value  of  the  inelastic  cross  section  at  high 
energies  (E  >  100  Mev).  Experimentally,  it  has  been  found  that  an  incident  proton 
flux  will  induce  the  emission  of  neutrons,  photons,  secondary  protons  and,  for  energies 
of  greater  than  300  Mev,  the  entire  spectrum  of  other  elementary  particles.  It  has 
been  shown  (Reference  84)  that  the  emission  of  particles,  other  than  the  first  three 
mentioned,  is  negligible  for  the  incident  spectra  of  interest  (solar  particle  events 
and  trapped  radiation  belts);  for  absorbed  dose  calculations,  they  can  be  neglected. 

In  the  special  case  of  a  hydrogen  absorber,  only  elastically  scattered  protons  can 
result  from  the  proton-proton  collision  for  incident  proton  energies  up  to  the  pion  pro¬ 
duction  threshold  near  300  Mev.  The  only  secondary  particles  emitted  in  this  case 
will  be  protons. 

Assume  a  flux  of  protons,  $  (E,x),  at  a  point  x  in  an  absorber.  Then  the  rate 
at  which  inelastic  collisions  are  JccurrJng  is  given  by 


r  max 

T?)  VE'X)'in(E)dE 
E  . 
min 


where 


(E,x)  *  incident  proton  flux  at  point  x 
P 

a.  *  inelastic  collision  cross  section 
in 

E  .  ,  E  ■  the  energy  ranges  of  the  proton  flux 
min  max 

The  particle  emission  processes  induced  by  these  collisions  are  energy  dependent. 
The  high-energy  protons,  which  are  defined  as  having  an  energy  much  larger  than  the 
binding  energy  of  the  individual  nucleons  of  the  struck  nucleus,  are  assumed  to  interact 
with  individual  nucleons  by  elastic  collisions.  This  intranuclear  cascade  process  has 
been  studied  extensively  by  computer  simulation  of  individual  collisions  (References  85^ 
86,  and  87),  the  Monte  Carlo  technique.  When  no  cascade  particle  remains  with 
sufficient  energy  to  escape  the  potential  field  of  the  nucleus,  the  intranuclear  cascade 
process  is  considered  terminated,  and  the  evaporation  model  is  used  to  further  predict 
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the  behavior  of  the  excited  nucleus  (References  88,  89,  and  90).  A  proton  with 
insufficient  energy  to  initiate  a  cascade  is  considered  to  excite  the  entire  nucleus, 
and  the  evaporation  theory  is  used  to  predict  the  decay  products.  These  theoretical 
methods  have  both  strengths  and  weaknesses,  and  detailed  comparisons  with  the 
available  experimental  data  are  available  in  the  cited  literature.  While  it  would  be 
most  desirable  to  use  only  experimental  results  for  the  analysis  of  the  induced  secondary 
radiations,  the  sparse  data  available  and  the  enormous  amount  of  information  required 
make  this  impossible  at  present.  Theoretical  predictions  must  be  utilized  to  fill  the 
gaps  in  the  present  experimental  data.  This  is  rather  unfortunate,  as  the  present 
knowledge  of  the  structure  of  the  nucleus  and  the  nature  of  nuclear  reactions  is  much 
less  precise  than  our  knowledge  of  the  proton-electron  collision  processes. 

The  following  data  must  be  obtained  to  quantitatively  describe  the  emission 
process: 

1)  The  inelastic  cross  section  of  the  absorber.  These  data  are  generally  avail¬ 
able  from  experiment  for  most  shielding  materials  of  interest.  Interpolation  on  both 
energy  and  atomic  weight  is  possible. 

2)  The  number  of  neutrons  or  protons  emitted  per  inelastic  collision  of  an 
incident  neutron  or  proton.  While  experimental  values  are  given  in  the  literature  for 
particular  materials  and  incident  particle  energies,  theoretical  calculations  provide 
the  bulk  of  available  data. 

3)  The  energy  distribution  of  the  emitted  particles.  Theoretical  calculations 
provide  a  direct  prediction  of  the  energy  distribution,  while  the  experimental  data 
generally  give  only  an  indirect  comparison  with  theory. 

4)  The  number  and  energy  distribution  of  gamma  rays  emitted  by  the  struck 
nucleus.  Here  there  is  an  abundance  of  experimental  data,  though  little  of  it  is  in  a 
form  usable  for  shielding  calculations.  As  the  gamma  rays  are  indicative  of  the 
energy  levels  of  the  nucleus,  they  have  been  used  for  the  analysis  of  nuclear  structure. 
The  experimental  results  of  Wakatsuki,  et  al.  (Reference  91)  are  in  a  convenient  form 
for  this  calculation,  and  his  results  have  been  used. 

The  proton  secondary  calculations  have  been  formulated  in  such  a  way  that 
these  basic  data,  obtained  from  either  theoretical  or  experimental  work,  can  be 
used  as  input  information.  The  results  reported  here  use  the  data  of  Wallace  and 
Sondhaus  (Reference  92),  which  are  essentially  based  on  the  Metropolis  calculations 
(References  85  and  86).  The  Oak  Ridge  group  has  recently  published  its  calculations 
in  a  form  that  is  directly  compatible  with  these  calculations  (Reference  93),  and 
either  source  may  be  used. 
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6.2  RADIATION  TRANSPORT 


Once  the  details  of  the  physical  interactions  between  the  incident  protons  and 
the  absorber's  atoms  are  specified,  either  by  theoretical  or  experimental  results,  it  is 
then,  in  principle,  possible  to  determine  the  radiations  emitted  by  the  absorber. 

While  idealized  absorber  geometries  and  simplified  physical  interactions  sometime 
permit  exact  solutions  for  the  radiation  transport,  more  realistic  problems  generally 
force  the  use  of  various  approximations.  It  is  desirable  in  shielding  calculations  to 
make  these  approximations  in  such  a  way  that  a  conservative  estimate  of  the  absorbed 
dose  is  made;  yet  the  results  must  be  as  accurate  as  possible  or  they  lose  their  useful¬ 
ness.  Fortunately  for  the  shielding  situations  currently  of  greatest  interest,  little 
accuracy  need  be  sacrificed. 

The  standard  transport  equations  have  been  used  by  workers  at  Oak  Ridge 
(Reference  84),  Lockheed  (Reference  94),  and  The  Boeing  Company  (Reference  95) 
to  describe  the  propagation  of  an  incident  proton  spectrum,  and  the  resulting  cascade 
and  evaporation  products  through  various  shield  configurations.  The  approach 
generally  used  in  the  solution  of  these  equations  (some  form  of  finite  difference  approxi¬ 
mation)  has  the  disadvantage  of  requiring  considerable  computer  time  to  achieve 
results  of  reasonable  accuracy.  By  formulating  the  transport  equations  as  a  series  of 
integral  equations,  it  was  found  possible  to  achieve  a  high  degree  of  accuracy  with  a 
minimum  of  computer  time. 

6. 2.  1  Protons 

6. 2.  1.1  Incident  Protons 

Consider  a  proton  spectrum  incident  normally  upon  an  absorber  of  finite  thick¬ 
ness  and  infinite  extent.  Assume  the  following  definitions: 

d$  the  proton  flux,  as  a  function  of 

-gjj-  “  $  (E,  x)  energy  and  position  in  the  absorber 

^  in  units  of  protons/cm  -Mev 

The  total  number  of  protons  is  then 

♦  W*  (*6(E,x)dE 
P  J  P 
o 

and  definitions  of  total  and  average  energy  of  the  proton  flux  follow  naturally.  The 
basic  problem  is  to  find  the  proton  flux  at  any  point  in  the  absorber 

♦p(E,x) 
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Recall  that  the  Incident  protons  travel  In  essentially  straight  lines  until  they 
are  stopped  or  make  an  inelastic  collision  with  a  nucleus  of  an  absorber  atom.  The 
energy  loss  per  unit  path  length  (-dE/dx)  can  then  be  used  to  calculate  the  total  energy 
loss  in  passing  through  the  absorber  by  equating  the  proton  path  length  and  the  thick¬ 
ness  of  the  absorber.  We  then  obtain  the  range-energy  relation  for  protons  in  terms 
of  the  initial  and  final  energy  of  the  proton  and  the  thickness  of  absorber  traversed. 


Then 


R(E.  )  -  x 

in 


d  R(E  ) 

out 


dE 


out 


d  R(E,n) 

~3T 


or 


dE 


dE. 


out  "in 

wimwn 

our  in 


S(E)- 


dE 

dx 


If  no  protons  are  stopped  in  the  absorber,  the  number  entering  must  equal  the  number 
emerging. 


d<}>  *  q>  (E.  ,  0) dE.  *  $  (E  .x)dE 
Yp  Tp  in  in  p  ouf  out 

where  E.  and  E  x  are  related  by  the  range-energy  relation, 
in  out  ' 

The  proton  spectrum  at  point  x  is  then  related  to  the  incident  spectrum  by 
$  (E  „x)«4>  (E.  ,  0)  S(Eln) 

V  ouf  V  OT;J-T 

v  out' 


The  proton  spectrum  at  any  point  in  the  absorber  can  be  obtained  in  this  manner. 
For  a  multilayer  absorber  of  different  material  types,  the  procedure  is  to  repeat  the 
calculation  as  the  proton  flux  passes  through  each  layer. 

For  thin  absorbers,  that  is,  x  much  less  than  the  mean  free  path  of  protons  in  the 
absorber  between  inelastic  collisions,  this  procedure  gives  proton  fluxes  of  adequate 
accuracy;  as  absorbers  of  increasing  thickness  are  considered,  it  becomes  necessary 
to  attenuate  the  incident  proton  flux  to  account  for  those  protons  removed  by 
inelastic  collisions.  When  this  is  done,  we  find 
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where 


|j(E* '  TT  "in.la.H, 


inelastic 


inelastic  cross  section  of  material 


By  using  the  relation 


to  obtain 


R(E)  -  R(E.n)  -  x 


dE  -  -  S(E)  dx 


the  variable  of  integration  is  changed  from  x  to  E, 

6. 2. 1.2  Secondary  Protons 

The  incident  proton  flux  generates  secondary  protons  in  the  absorber  through 
inelastic  collisions  with  absorber  nuclei.  While  these  secondary  protons  may  be 
emitted  with  angular  distributions,  ranging  from  isotropic  to  highly  anisotropic,  in 
what  follows  they  will  be  considered  as  traveling  in  the  direction  of  the  incident 
proton.  This  approximation  will  be  discussed  later.  For  an  absorber  of  thickness  X, 
the  generation  rate  of  secondary  protons  of  energy  E1  at  point  x  in  the  absorber  can 
be  given  as  follows: 

E 

n  r max 

-s2  j  vE'x)vE'M)FP<E'E,)dE  dx 

^  J  E' 

where  all  terms  are  defined  as  follows: 

23 

N  *  Avogadro's  number  “6. 025  (10  )  atoms/mole 
o 

A  *  atomic  weight  of  material  M,  grams/mole 

-2  -1 

f  (E,x)  =  proton  flux  at  energy  E  and  point  x  in  units  of  proton  cm  Mev 

P 
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°  (E,M) 
PP 


yield  cross  section  for  the  production  of  protons  by  the  inelastic 
collision  of  protons  of  energy  E  with  a  nucleus  of  element  M 


a  (E,M)-N  (E,M)<7  (E,M) 
pp  p  in 

N  (E,M)  ■  number  of  protons  per  inelastic  collision  of  protons  with  energy  E 
p  on  nuclei  of  material  M 


o. 

in 


(E/M) 


inelastic  collision  cross  section  in  units  of  cm  /atom 


FpM') 


spectral  distribution  of  protons  produced  by  protons  of  energy 
E  in  units  of  proton/Mev  proton 


Set 


E  ■  maximum  energy  of  the  incident  proton  flux 


E'  «  energy  of  secondary  proton 
E 


N 


PP 


fmax 

♦r* 


x)  a  (E,M)  F  (E,E')  dE 
PP  P 


Now  the  flux  of  secondary  protons  at  energy  E  emerging  from  the  shield  of 
thickness  Xi,  * 


W-fw^fy 

o  ,p 

where  E  and  E*  are  related  by 

sp  7 

R(E')  -  R(E$p)  +  (X-x) 

the  range-enei^y  relation  for  the  absorber*  While  this  flux  of  secondary  protons  in 
turn  is  further  attenuated  by  inelastic  collisions  and  produces  third-generation  par¬ 
ticles,  the  process  will  be  terminated  at  this  point.  The  results  indicate  that  for  the 
proton  spectra  and  shield  materials  of  interest  in  space  vehicle  shielding,  no 
appreciable  error  is  made  in  terminating  the  cascade  at  the  second  order. 


6. 2. 2  Neutrons 


The  incident  proton  flux  also  generates  neutrons  in  the  absorber.  These 
neutrons  are  moderated  in  energy  by  elastic  collisions  with  light  nuclei,  hydrogen 
being  the  most  effective.  Inelastic  collisions,  which  effectively  remove  the  neMtron 
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from  the  particle  flux,  are  here  considered  to  be  the  only  attenuation  mechanism.  In 
those  cases  where  more  precise  neutron  flux  solutions  are  desired,  this  calculation  can 
be  used  to  obtain  the  neutron  source  terms.  The  neutron  is  assumed  to  travel  in  the 
direction  of  the  incident  proton. 


For  a  single  absorber,  the  source  distribution  of  neutrons  of  energy  E  can  be 
written  as  follows* 


Now  let 


N 

o 

X 


E 

-  max 


lVE'X>  apn(E,M)  Fn  (E'E‘)dE 
E‘ 


dx 


P 

pn 


(E',x) 


N 

o 

X 


E 

N  max 

U  (E,x)  a  (E, M)  F  (E,  E')  dE 
J  P  Pn  n 

E' 


The  symbols  are  defined  as  follows: 


23, 


N  »  Avogadro's  number  *  6.  025  (10  )  atoms/mole 
o 

A  «  atomic  weight  of  material  M,  grams/mole 

-2  -1 

<^(E,x)  *  proton  flux  .at  energy  E  and  point  x  in  units  of  protons  cm  Mev 

o  (E,M)  *  yield  cross  section  for  the  production  of  neutrons  by  the  inelastic 
collision  of  protons  of  energy  E  on  a  nucleus  of  element  M 

a  (E,M)  -  N  (E,M)  a.  (E,M) 
pn  n  in 

N  *  number  of  neutrons  per  inelastic  collision  of  protons  of  energy  E 
with  atoms  of  material  M 


*  inelastic  collision  cross  section  in  units  of  cm  /atom 

F  (E,  E')  *  spectral  distribution  of  neutrons  produced  by  protons  of  energy  E 
(In  units  of  neutrons/Mev  proton) 


max 


maximum  energy  of  incident  proton  flux  (Mev) 


Integrating  the  neutron  generation  function  P  (E,x)  over  the  shield  thickness 
X  gives  the  total  number  of  neutrons  of  energy  E  generated  in  the  absorber. 
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X 


<t>  (E'M  P  (E',x)  dx 
n  J  pn 

o 

if  it  it  assumed  that  all  the  neutrons  generated  travel  in  the  direction  of  the 
incident  protons,  and  that  those  neutrons  that  undergo  inelastic  collisions  are  effec¬ 
tively  removed  from  the  neutron  beam,  the  number  of  neutrons  both  generated  and 
penetrating  the  absorber  is 


X 

<t>n  (E'/X)  m\ Ppn^,/X)  exp[-pn(E,,M)(X-x)]  dx 


where  u  (E',M)  is  the  attenuation  coefficient  for  neutrons  of  energy  E1  and  material  M 
In  cm /gm. 

6. 2. 3  Gammas 


The  inelastic  collisions  of  the  incident  protons  with  absorber  nuclei  lead  to 
gamma  emission  by  the  excited  nuclei.  These  gamma  rays  are  quite  penetrating, 
and  their  contribution  to  the  absorbed  dose  must  be  evaluated. 


The  rate  at  which  the  incident  protons  are  suffering  inelastic  collisions  in  the 
absorber  may  be  expressed  as  follows: 


max 


<|>  (E,x)  a.  (E,M)  dE 
p  in 


where  the  symbols  are  as  defined  previously. 


Integrating  over  the  absorber  thickness  and  including  a  transmission  factor  for 
Hie  attenuation  of  the  gammas  gives  the  following  expression  for  the  flux  of  gamma 
rays  at  the  surface  of  the  absorber,  assuming  all  the  gammas  are  emitted  in  the  direc¬ 
tion  of  the  Incident  proton. 


♦y<P> 


p%  r 


max 


*p(E,x)  *Jn(E,M)  dE 


®xpE-Hy(E',  M)(X-x)]  dx 


126 


where 


<j>^(E')  *  photon  flux  at  energy  E1  in  photons/Mev 

Fy(E',M)  *  spectral  density  of  photons  of  energy  E'  from  material  M 

N  *  Avogadro’s  number 
o 

A  *  atomic  weight  of  the  material 

a.n(E,M)  *  inelastic  cross  section  of  nuclei  of  type  M  at  incident  proton  energy  E 
u  (E',M)  *  attenuation  coefficient  for  gammas  of  energy  E'  in  material  M 

r 

X  m  total  thickness  of  the  absorber 


6.2.4  Geometrical  Considerations 


By  the  solution  of  the  radiation  transport  equations,  the  radiation  flux  behind 
a  multilayer  slab  shield  can  be  found.  To  relate  the  one-dimensional  shielding 
problem  considered  by  the  transport  equations  to  a  realistic  three-dimensional  space 
vehicle,  the  angular  distributions  of  the  primary  and  secondary  radiations  must  be 
considered.  The  incident  proton  flux  is  defined  by  an  isotropic  flux  intensity  J  . 
This  is  the  number  of  protons  per  second  that  pass  through  a  hypothetical  sphere  of 
1  cm^  cross-sectional  area.  The  number  that  enters  the  sphere  from  a  given  solid 
angle,  AQ,  is  then 


AN 


4  “ 


It  can  be  shown  that  the  average  path  length  of  a  particle  striking  the  hypo¬ 
thetical  sphere  is  such  that  if  a  monodirectional  flux  and  an  isotropic  flux  of  the  same 
intensity  (particles/cni  )  are  incident  on  an  absorber,  the  resulting  doses  are  equal. 

If  the  incident  particles  travel  essentially  in  straight  lines  through  the  absorber,  the 
three-dimensional  vehicle  analysis  can  be  reduced  to  a  series  of  one-dimensional  slab 
shielding  problems. 

When  the  generation  of  secondary  radiation  is  considered,  their  angular  distri¬ 
butions  pose  a  complex  problem.  The  basic  difficulty  is  not  only,  as  one  might 
suppose,  finding  an  adequate  representation  for  the  respective  angular  distributions 
as  a  function  of  particle  energies,  but  also  finding  the  angular  distribution  of  the  pro¬ 
tons  generating  secondaries  in  the  absorber.  To  accomplish  this  would  require  a 
three-dimensional  treatment  of  the  entire  vehicle,  a  prohibitively  lengthy  analysis 
for  any  realistic  space  vehicle.  To  avoid  this,  the  angular  distribution  of  each 
radiation  type  will  be  analyzed  and  a  simplified  correction  term  will  be  selected. 
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The  angular  distribution  of  the  secondary  protons  depends  on  the  ratio  of  the 
emitted  secondary  to  the  incident  primary  proton.  Collision  kinematics  require  that 
secondaries  of  relatively  high  energy  be  emitted  with  a  distribution  highly  peaked  in 
the  direction  of  the  incident  particle.  As  lower-energy  particles  are  considered,  the 
distribution  gradually  broadens  to  cover  the  entire  forward  hemisphere.  Evaporation 
particles  are  emitted  isotropically. 

As  only  the  high-energy  particles  have  appreciable  range,  and  as  the  law-energy 
particles  deposit  all  their  energy  locally,  assume  that  all  secondary  protons  are 
transmitted  in  the  forward  direction  —  the  straight-ahead  approximation. 

The  angular  distribution  of  the  emitted  neutrons  is  similar  to  that  of  the  secondary 
protons.  Because  of  the  greater  penetrating  power  of  the  neutrons,  it  would  be 
incorrect  to  assume  local  absorption  of  the  low-energy  neutrons.  A  rigorous  treatment 
of  neutron  transport  would  involve  a  complete  consideration  of  the  entire  space  vehicle. 
As  this  would  involve  a  very  time-consuming  analysis,  for  the  purposes  of  this  program 
it  will  be  assumed  that  all  neutrons  are  transmitted  in  the  direction  of  the  incident 
proton.  The  gamma  rays  are  also  treated  by  the  straight-ahead  approximation.  These 
approximations  insure  that  the  secondary  dose  calculation  will  give  an  upper  limit  in 
nearly  all  situations.  For  a  strongly  anisotropic  incident  flux,  special  configurations 
of  shielding  can  invalidate  this  approximation  and  these  situations  should  be  treated 
with  caution. 

6.3  DOSE  CALCULATIONS 
6.3.1  Energy  Deposition 

When  the  spectra  of  incident  protons  and  their  induced  radiations  have  been 
found  at  any  point  in  the  astronaut's  body  or  equipment,  the  absorbed  dose  can  be 
found.  The  unit  of  absorbed  dose  used  Is  the  rad,  defined  as  100  ergs  of  absorbed 
energy  per  gram  of  material.  The  crucial  element  in  the  definition  is  that  the  energy 
actually  be  absorbed  in  the  material.  To  calculate  the  absorbed  dose  in  an  accurate 
manner,  it  is  necessary  to  know  not  only  the  energy  loss  to  the  absorber  through  pri¬ 
mary  collisions  of  the  incident  radiation,  but  also  the  secondary  and  higher-order 
radiations  created  by  the  primary  energy  loss  must  be  followed  until  their  energy  is 
low  enough  to  ensure  that  their  remaining  energy  will  be  deposited  locally,  or  that 
they  have  left  ihe  area  of  interest.  The  general  result  of  this  requirement  is  that  the 
absorbed  dose  at  a  particular  point  in  an  absorber  depends  not  only  on  the  primary 
radiation  spectrum  at  that  point,  but  also  on  the  secondary  radiations  generated  in 
the  entire  absorber  that  penetrate  to  the  point  of  interest. 

The  incident  proton  flux  deposits  energy  In  on  absorber  by  ionization  of  the 
absorber  electrons.  The  energy  lost  per  gm/cm2  is  given  by  the  Beflhe-Bloch 
equation  and  can  be  used  to  find  the  energy  lost  in  the  following  way| 
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These  electrons  have  then  a  very  limited  range  for  the  Incident  proton  energies 
of  Interest  (<1000  Mev),  and  their  total  Ionization  energy  loss  is  quite  localized^  The 
secondary  protons  may  be  treated  in  the  same  manner. 


The  neutrons  generated  by  the  incident  protons  present  a  more  complex  problem. 
As  the  neutron  is  not  a  charged  particle,  it  deposits  energy  only  by  collisions  with  the 
absorber's  nuclei.  These  collisions  in  turn  result  in  a  transfer  of  energy  to  the  struck 
nucleus,  either  by  absorption  of  the  incident  neutron,  or  by  elastic  or  inelastic 
scattering.  In  any  case,  the  initial  energy  of  the  neutron  can  be  dissipated  in  a 
variety  of  ways,  ranging  from  emission  of  a  y-r ay,  neutron,  proton,  or  heavier  par¬ 
ticle  by  the  excited  nucleus  to  the  ionization  in  the  absorber  resulting  from  the  kinetic 
energy  impasted  to  the  nucleus  in  the  collision.  An  exact  calculation  of  the  absorbed 
dose  would  follow  these  tertiary  particles  through  their  respective  histories  as  energy 
depositors.  As  this  would  take  a  great  deal  of  ttme,  an  approximate  dose  conversion 
factor  developed  by  W.  A.  Gibson  (Reference  96)  is  used  which  gives  an  upper-limit 
estimation  of  the  absorbed  dose. 


For  space  vehicle  shielding  analysis,  where  thin  shields  and  low  secondary  doses 
are  the  rule,  this  method  provides  sufficient  accuracy.  The  dose  at  a  point  x  Is  then 
expressed  as  follows: 

E 

r  max 

D to -\  9  (E,x)D  (E)dE 
J  n  n 

^mln 


where 


$n(E,x)  *  neutron  flux  at  point  x  In  neutron/Mev  cmi 


E  m  maximum  neutron  energy  considered 
max 


E  .  *  minimum  neutron  energy  considered 
min 
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D 


(E)  *  neutron  flux  to  dose  conversion  factor  in  rads/neutron/cm 
n 


2 


D(x)  *  dose  in  rads  at  point  x 


The  gamma  rays  present  a  problem  similar  to,  though  less  complex,  than  the 
neutrons.  The  primary  absorption  mechanisms  of  high-energy  gammas,  Compton 
scattering  and  pair  production,  result  in  energetic  electrons  and  photons  of  degraded 
energies.  The  range  of  the  electrons  is  small,  and  their  energy  can  be  considered  as 
deposited  locally.  The  secondary  gammas,  however,  have  mean  absorption  ranges 
that  can  be  significant  in  comparison  to  the  human  body.  Due  to  their  relative  lack 
of  importance  in  comparison  to  the  other  radiations,  a  simple  upper-limit  dose  con¬ 
version  factor  will  be  used.  The  expression  for  the  gamma  dose  then  has  the  same 
form  as  the  other  dose  expressions, 

E 

r  max 

E>rM-3  ♦r(E,x)Dj,(E)  dE 

E  . 

min 


6.3.2  Linear  Energy  Transfer  (LET)  Spectrum 

The  response  of  a  biological  system  to  radiation  has  been  shown  to  depend  not 
only  on  the  total  energy  deposited  —  the  absorbed  dose  — but  also  on  the  rate  at  which 
the  radiations  deposit  energy  to  the  tissue  they  traverse,  the  LET,  or  dE/dx.  A 
charged  particle  has,  for  any  energy,  a  given  dE/dx,  available  either  from  the 
Bethe-Bloch  equation,  or  from  experimental  results.  A  monoenergetic  beam  of  pro¬ 
tons  incident  upon  a  tissue  slab  would  have  a  particular  energy  and  dE/dx  at  each 
point  in  the  tissue,  neglecting  beam  straggling.  An  incident  spectrum  of  protons 
would,  of  course,  produce  a  spectrum  of  dE/dx  values  at  any  point  in  the  absorber. 

From  the  previously  discussed  proton  transport  relations,  the  proton  flux  at  a 
point  x  in  the  absorber  can  be  obtained.  From  this  a  transformation  of  variables  is 
made: 


d*p  *  *p  x) dE  *  f p(y» x)  dy 


where 

dE 


dy 


d 


dE 

dx 
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now 


The  relation  between  E  and  dE/dx  Is  given  by  the  Bethe-Bloch  equation.  An 
empirical  expression  can  be  used  to  fit  the  dE/dx  versus  E  plot,  thus  permitting  rapid 
evaluation  of  the  derivative,  and  in  addition,  extending  the  range  of  valid  LET  spec¬ 
trum  determinations  below  the  limits  of  the  Bethe-Bloch  equation.  The  secondary 
proton  spectrum  can  be  evaluated  in  the  same  way.  One  precaution  must  be  observed. 
For  the  straightforward  transformation  from  $  (E,x)  to  ^  (dE/dx,  x)  to  be  valid  one 
must  operate  in  a  region  where  E  is  a  single-valued  function  of  dE/dx.  For  tissue, 
the  region  of  0.  1  to  1000  Mev  satisfies  this  requirement,  and  is  also  the  region  of 
primary  interest  In  dose  calculations. 

Once  the  number  LET  expression  (])  (dE/dx,  x)  has  been  obtained,  it  can  be 
used  to  determine  some  quantities  of  interest.  For  example,  the  number  of  protons 
between  any  two  values  of  dE/dx  is  given  by 


dE 


the  energy  deposited  by  protons  having  dE/dx  values  in  the  range  from  (dE/dx).  to 
(dE/dx)^  is  then  * 


dE 


2  _i 

in  our  calculations  the  following  units  are  t/sed:  dE/dx  is  in  Mev  cm  gm  , 

5  (dE/dx,  x)  is  In  (protons/cm  )/(Mev  cm  gm  )  and  E  is  in  Mev.  To  convert 
from  the  density  independent  form  dE/dx  to  the  common  radiobiological  units  of  kev 
per  micron,  a  density  for  tissue,  or  the  absorber  in  question,  must  be  assumed.  One 
finds  that  the  conversion  factor  K  from  Mev  cm V gm  to  kev/micron  is 

k*  a  1  pt 


kev/micron 

Mev  cm  /gm 


131 


3 

where  p ^  is  the  density  of  tissue  in  gm/cm  ,  usually  taken  to  be  unity. 
6. 3.  3  Semiconductor  Damage 


In  addition  to  the  effect  of  absorbed  energy  on  the  human  body,  semiconductor 
devices  have  also  been  found  sensitive  to  radiation  damgge.  The  defects  produced 
by  an  incident  radiation  are  related  to  the  operating  parameters  of  various  solid-state 
electronic  components  (Reference  97).  As  the  defect  density,  defects  per  cm  ,  is  an 
energy-dependent  function  of  the  incident  protons,  it  can  be  calculated  in  a  manner 
similar  to  the  absorbed  doses 


DD(x) 


E 

s  max 


)  ♦  (E,x)  D(E)  dE 


'min 


where 


(}>p(E,x)  *  proton  spectrum  at  point  x  in  proton^cm- Mev) 

D(E) »  (dislocations/cm^)  per(proton/cm^) 

3 

DD(x)  »  defect  density  at  x  in  defects/cm 

The  function  D(E)  has  been  deduced  from  a  variety  of  experimental  results  on 
silicon,  shown  in  Figure  44.  Similar  curves  for  germanium  will  be  available  when  the 
experimental  data  presently  available  are  reduced. 

6.4  RESULTS 


6.4.1  Program  Descriptions 

Two  computer  programs  have  been  developed  to  evaluate  the  radiation  transport 
equations  and  to  perform  the  dose  calculations  described  for  incident  protons.  The 
computational  techniques  used  to  solve  these  equations  and  detailed  operating  instruc¬ 
tions  are  given  in  Volume  II  of  this  report.  One  program,  the  proton  primary-dose 
code,  was  developed  for  the  evaluation  of  the  absorbed  dose  resulting  from  the  direct 
ionization  energy  loss  of  penetrating  protons.  The  number  LET  and  energy  LET  spectra, 
as  well  as  the  local  dose,  are  also  evaluated.  The  basic  approximation  made  in  this 
code  is  the  neglect  of  proton  inelastic  collisions.  It  will  be  shown  that  this  approxi¬ 
mation  is  of  adequate  accuracy  for  the  thin  shields  of  present  space  vehicles  and  most 
incident  proton  fluxes.  Only  for  the  very  soft  spectra  are  the  secondary  radiations 
important  for  thin  shields.  The  proton  secondary  code  evaluates  the  contribution  made 
by  the  secondary  radiation  to  the  total  dose.  The  number  LET  spectra  of  the  primary 
and  secondary  protons  are  also  calculated.  The  proton  primary  code  has  the  advantage 
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FIGURE  44.  Density  of  dislocation  formation  in  silicon  due  to  monoenergetic  protons 


of  great  speed,  allowing  many  combinations  of  shielding  and  incident  proton  spectra 
to  be  evaluated  in  a  minimum  time.  The  run  time  T  in  minutes  per  sector  for  each 
shielding  combination  or  sector  is  approximately 

T-Q.03  N Tjg 

where 


N  *  number  of  layers  of  material  in  the  sector 
M  *  number  of  output  energy  points,  usually  50 


The  more  accurate  proton  secondary  code  has  a  rtin  time  of 
T-0.4N| 


6.4.2  Sample  Calculations 

To  demonstrate  the  geneial  features  of  the  calculations,  some  representative 
results  are  given.  First,  in  Figures  45  through  48  are  shown  curves  of  dose  versus 
aluminum  thickness  for  a  trapped  proton  spectrum  and  three  solar  particle  event 
spectra.  The  relative  importance  of  the  prtmary  and  secondary  dose  components  is 
heavily  dependent  on  the  spectral  form  of  the  incident  proton  flux.  The  regton  of 
validity  of  the  proton  primary  code  is  a  function  of  the  "hardness"  of  the  incident  pro¬ 
ton  spectrum.  It  can  be  seen  that  with  the  exception  of  the  40-Mv  solar  particle  event 
spectrum,  the  proton  primary  code  will  give  reasonable  dose  estimates  for  even  very 
thick  shields.  The  two  outstanding  features  of  these  dose  curves  are  the  very  rapid 
decrease  in  dose  in  the  first  few  gm/cm^  of  absorber,  and  the  gradual  transition  to 
a  very  slaw  attenuation  with  thickness  region  as  the  low-energy  protons  are  removed 
and  the  more  penetrating  secondary  radiations  begin  to  be  felt. 


The  number-LET  spectra  plotted  in  Figures  49  through  51  for  representative 
solar  particle  events  show  that  a  few  gm/cnv  of  shielding  quickly  harden  the  LET 
spectra.  The  energy-LET  spectra  given  in  Figure  52  show  that  even  for  the  softest 
incident  proton  spectrum,  a  few  gm/citP  of  shielding  shifts  the  peak  energy  deposition 
region  to  low  LET  values. 


The  differential  energy  spectrum  of  protons  from  a  100-Mv  particle  event  passing 
through  various  thicknesses  of  aluminum  are  shown  in  Figure  53.  The  proton  spectrum 
quickly  assumes  a  constant  spectral  shape  which  gradually  hardens  as  increasing 
thicknesses  are  reached. 
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FIGURE  46.  40-Mv  tolar  avant  spoctrum  —  protons  Incident 
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Protons/cm  per  kev/micron 


FIGURE  49,  Number-LET  spectrum  for  40-Mv  solar  event 


kev/micron 


FIGURE  50,  Number- LET  spectrum  for  100-Mv  solar  event 


Proton*/cm  per  kev/micron 
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6,4.3  Comparison  with  Oak  Ridge  Calculations 

In  Figure  54,  a  comparison  between  the  results  of  Alsmiller  and  Murphy 
(Reference  84)  and  the  present  proton  secondary  calculation  is  given.  Although  the 
calculations  were  made  as  comparable  as  possible,  certain  differences  remain.  The 
results  of  Alsmiller  are  given  for  particles  of  energy  32  Mev  or  greater,  and  the 
incident  proton  spectrum  has  a  finite  cutoff  at  1000  Mev.  In  the  present  calculations, 
the  lower  energy  limit  is  30  Mev,  and  the  incident  spectrum  is  extended  past  1000 
Mev  by  extrapolation.  The  nuclear  input  data  of  Alsmiller  was  obtained  from  cosmic 
ray  measurements,  while  our  results  were  based  on  the  data  of  Wallace  and  Sondhaus. 
Because  of  these  differences  in  the  calculations,  it  is  felt  that  the  results  are  in 
satisfactory  agreement. 
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Aluminum  thickness  (gm/cm^) 

FIGURE  54.  Comparison  of  the  results  of  Alsmiller  with  this  study 


7.0  ANGULAR  DISTRIBUTION  OF  INCIDENT  RADIATION 
FLUXES  AND  THE  EFFECT  OF  VEHICLE  ORIENTATION 


In  previous  space-radiation  shielding  analyses  (References  42  and  98),  the 
incident  flux  was  assumed  to  be  isotropic.  This  study  was  made  to  establish  the  per¬ 
turbation  on  the  isotropic  flux  results  when  the  incident  flux  is  nonisotropic.  All 
cases  studied  assumed  that  the  specified  incident  flux  is  symmetric  about  an  axisTP 
and  is  dependent  only  upon  the  angle  6,  measured  abouto,  and  not  the  angle  t|>,  as 
shown  in  Figure  55. 


FIGURE  55.  Definition  of  angles 

Three  types  of  particle  deflections  were  studied:  (1)  transmission  where  an 
incident  particle  is  degraded  in  energy  but  retains  its  original  direction,  (2)  trans¬ 
mission  where  the  incident  particle  is  degraded  and  the  distribution  in  final  direction 
is  described  by  f(0),  and  (3)  transmission  of  secondary  particles  of  definite  directional 
distribution. 

In  the  first  case,  the  vehicle  was  characterized  by  a  fixed  coordinate  system  as 
shown  in  Figure  56.  The  transmission  where  the  initial  direction  was  maintained  after 
a  collision  was  typical  of  proton  penetration. 

Three  positions  of  B,  the  axis  of  symmetry  for  the  incident  flux,  were  studied 
(0,0,  1),  (0,707,  0,  0,707), and  (1,0,0)  (see  Figure  57),  and  three  differential  angular 
distributions  of  the  form 

f(0)  dO  *  A  sin  (0  -  0  )  d  for  0  >  0 

o  —  o 

*  0  for  0  <  0 

o 

where  A  is  a  normalization  constant  which  allows  the  integral  of  f(0)  about  4w 
steradians  to  be  equal  to  unity.  Values  of  30,  45,  and  60  degrees  were  used  for  0q. 

The  Boelr^  automatic  sector  analysis  program  (Reference  42)  was  used  to 
generate  the  sector  data  for  the  test  cylinder,  using  146  sectors.  Each  sector  was 
analyzed  to  find  the  slant  thickness  in  the  sector,  the  angle  0  between  the  incident 
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FIGURE  56.  Coordinate  system  for  angular  distribution  calculations 
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particle  path  and  B,  the  angle  ^  *=(ir/2j-  u  between  the  normal  to  the  sector  N 
and  B,  and  the  angle  y  between  tRe  incident  particle  path  and  the  normal  to  the 
sector  N.  The  angles  k  and  y  will  be  used  for  the  electron  and  bremsstrahlung 
analysis  and  are  furtherraescribed  in  Figure  58.  The  receiver  point  for  this  study 
was  defined  as  the  center  of  the  cylinder. 


The  sector  data  was  used  as  input  for  the  proton  code  and  the  total  dose 
transmission  factor  was  computed  from  the  relation 


D  * 
P 


I 

y~  u.  f.  (9)  t. 

* —  ii  i 

i*=  1 


where  u,  is  the  fractional  solid  angle  subtended  by  sector  i,  f.  (0)  is  the  fraction  of 
incident  radiation  entering  sector  i  traveling  towards  the  dose  point,  and  T.  is  the 
transmission  of  the  incident  flux  through  sector  i.  It  is  assumed  in  this  analysis  that 
the  angular  distribution  of  the  incident  radiation  is  independent  of  the  incident 
energy  spectra. 

A  cylinder  is  drawn  to  represent  the  vehicle,  X is  the  reference  axis  (in  this 
case  possibly  the  spin  axis),  and  a  point  is  located  by  angles  w  measured  from  the^T 
and  an  azimuthal  angle  a.  If  the  vehicle  is  oriented  wjth  respect  to  the  axis  of 
symmetry  of  the  radiation  flux  so  that  u  and  a  define  B  relative  toX  the  number  of 
particles  N(AO)  incident  on  a  surface  element°dA  and  penetrating  to  the  origin  of 
the  coordinate  system  can  be  expressed  as: 

H(KQ)  *  N(0)  T(A) 


where 


cos  0  “  X*B 

N(0)  is  the  number  of  particles  PI  l"B|  incident  at  angle  0  on  area  dA  and  T(A)  is 
the  transmission  of  proton  through  the  shield  thickness  in  area  A. 

The  example  used  in  this  analysis  was  a  simple  orbiting  laboratory  in  the  shape 
of  a  right  circular  cylinder,  20  feet  long  and  130  inches  in  diameter.  Eight  wall 
thicknesses  studied  are  shown  in  Table  12. 
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FIGURE  58.  Geometry  for  cose  2  and  case  3  examples  of  nonisotropic 
angular  distribution 


Table  12 


SHIELD  CONFIGURATION  FOR  STUDY  OF  PROTON 
ANGULAR  DISTRIBUTION 


Case 

Side  Thickness 
(gm/cm^) 

End  Thickness 
(gm/cm^) 

1 

1 

1 

2 

2 

2 

3 

5 

5 

4 

0 

1 

5 

0 

2 

6 

1 

2 

7 

1 

5 

8 

2 

5 

The  results  of  this  study  are  shown  in  Table  13.  The  major  effect  of  the  non- 
i  so  tropic  distribution  was  the  redistribution  of  the  flux  through  sectors  of  greatest 
thickness  resulting  in  a  lower  dose  transmission  factor.  A  sine  distribution  produces 
a  dose  transmission  factor  20  percent  lower  than  the  isotropic  flux,  and  the  change 
in  cutoff  angle  from  30  to  60  degrees  had  a  very  small  effect.  For  the  proton  trans¬ 
mission,  these  results  indicate  that  a  nonisotropic  flux  distribution  of  the  sine  type 
will  produce  lower  transmission  factors  and  will  provide  a  safety  factor  when  isotropic 
flux  distributions  are  used  in  the  design  of  space  radiation  shields. 

Case  2  considers  the  electron  penetration  which  is  assumed  to  diffuse  through 
the  shield.  The  number  of  electrons/cm^  N.  ( £  )  crossing  the  surface  of  a  sector  is 
described  as  1 

2*  IT 

Nj  (  r  j"  F(0)  cos  |  sin  ^  d^  dB 

B-o 

For  an  isotropic  distribution,  F(0)  becomes  unity,  and  N.(^  )  becomes  $  /4  as  derived 
by  Evans  (Reference  64),  where  is  the  free-space  flux. 

The  fraction  of  the  free-space  flux  which  crosses  a  surface  can  be  computed  by 
integrating  the  fraction  entering  per  solid  angle  over  the  2*  solid  angle  above  the 
surface.  An  analytical  solution  could  not  be  obtained,  so  -iba'automatic  sector 
analysis  was  used  to  divide  the  2w  space  into  50  sectors  of  equal  solid  angle.  Using 
a  given  angle  ^  between  the  Taxis  and  the  normal  to  the  surface,  the  angle  6 
between  the  Incident  path  through  a  sector  and  Twas  computed.  The  fraction  N^(  j?  ) 
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penetrating  the  surface  can  be  expressed  as 


50 


F(0)  AQ.  cos 
i  Jo  l -  |  | 

ri 


where  F(0.)  is  the  fraction  of  the  flux  entering  through  solid  angles  AQ..  When  the 
angle  jr  is  zero,  the  expression  could  be  solved  analytically.  This  cise  was  used  as 
a  check  on  the  numerical  integration.  The  accuracy  was  better  than  10  percent. 
Results  of  the  analysis  of  N.(  j  ^  are  shown  in  Table  14.  Again  the  effect  of  a  non¬ 
isotropic  flux  was  to  reduce  the'' number  of  particles  transmitted,  since  in  the  dose 
expression 


D 

e 


I 


-L  T,  ♦,  N,  (  tj  F,W 

i-i 


N;(  f0)  will  be  less  than  the  isotropic  case  and  all  other  factors  are  unaffected  by  the 
angular,  distribution.  T.  is  the  transmission  through  sector  i,  (j>.  the  incident  free-space 
flux,  and  F.(y)  the  fraction  of  penetrating  particles  which  will1  be  intercepted  by  the 
receiver. 

The  transmission  programs  contain  at  least  two  dummy  variables  to  describe  the 
angular  behavior.  In  the  case  of  protons,  the  function  f.(0)  to  describe  the  fraction  of 
incident  flux  incident  on  sector  i  is  present.  For  brerrisstrahiung,  in  addition  to  N.(  f  ) 
and  F.(y),  expressions  describing  the  thick-slab  angular  distribution  as  a  function  of 
incident  and  exit  angle  can  be  employed.  At  the  time  of  this  writing,  there  are 
insufficient  data  to  employ  these  improvements;  however,  these  methods  can  be  used 
to  evaluate  hypothetical  situations. 

In  Case  3,  where  the  direction  of  the  secondary  is  related  to  the  direction  of 
die  incoming  particles,  the  problem  is  a  combination  of  Case  1  and  Case  2.  For 
bremsstrahlung  production,  each  sector  should  be  analyzed  to  determine  the  fraction  of 
bremsstrahlung  produced  which  will  be  directed  towards  the  receiver.  The  fraction  of 
photons  FBj(y)  is  expressed  as  (see  Figure  58) 


FB.(y)-  ^ G j (a)  G2(0)  sin  \  d£  dO 
A 

where  G.(a)  is  the  probability  to  scatter  toward  the  receiver,  and  G?(fl)  is  the 
probability  to  be  incident  at  an  angle  0  with  %  the  axis  of  symmetry?  The  dose  D 

B 
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Table  14 


FRACTION  OF  FREE  SPACE  FLUX  CROSSING  A  UNIT  SURFACE 


Angle  Between 

Normal  and 

Axis  of  Symmetry 

(degrees) 

Cutoff  Angle 
for  Sine 
Distribution 

(degrees) 

Fraction  of  Free 
Space  Flux  Crossing 
Unit  Surface 

0 

0 

0.21 

0 

30 

0.  13 

30 

30 

0.23 

45 

30 

0.20 

60 

30 

0.17 

90 

30 

0. 13 

0 

45 

0. 13 

30 

45 

0.27 

45 

45 

0.23 

60 

45 

0. 18 

90 

45 

0. 13 

0 

60 

0.048 

30 

60 

0.  14 

60 

60 

0.08 

90 

60 

0.048 

Isotropic  Flux 

0.25 
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will  be. 


v5L»ei  N(  L)ti  FBiw 


i-1 


where 

*  .  *  incident  electron  flux 

Tei 

N.(  *=  fraction  of  the  flux  crossing  the  surface 

T.  *  dose  transmission  factor 

i 

FB;(y)  *  fraction  directed  toward  the  receiver 

Provision  has  been  made  to  incorporate  G^(a)  and  FB.(y)  in  the  calculation  of  T.  (see 
Volume  II).  No  sample  calculations  were  made  in  this  case  since  G^(a)  could  not  be 
evaluated. 

The  problem  of  scatter  in  a  vehicle  was  briefly  examined;  however,  this  problem 
was  not  serious  enough  in  conventional  spacecraft  to  warrant  additional  code  develop** 
ment.  It  is  suggested  Jhat  for  protons,  no  correction  be  made  since  all  scattering  is 
predominantly  forward.  For  electrons,  the  backscattering  increases  with  atomic  num¬ 
ber  and  with  electron  energy.  Fortunately,  the  electrons  are  usually  almost  totally 
absorbed  in  the  vehicle  structure  (i.  e. ,  2  gm/cm^)  and  the  penetrating  particles  are 
degraded  to  an  extent  where  the  backscattered  fraction  is  a  negligible  quantity.  If 
a  very  thin  structure  is  encountered,  it  is  recommended  that  the  backscatter  correction 
factor  be  expressed  as 

J 

B  Ra* 

i-i 


where  B  is  the  enhancement  of  the  penetrating  flux  by  backscattering,  R  is  the  back¬ 
scatter  coefficient,  and  a  the  average  fraction  of  the  interior  surface  which  can  be 
viewed  from  a  point  on  the  interior  surface.  For  a  complicated  vehicle,  the  interior 
may  be  divided  into  several  areas,  such  as  the  crew  member,  the  visible  walls,  and 
the  subsystems.  For  Icw-Z  materials,  such  as  man  and  aluminum  structure,  the 
reflection  coefficient  will  be  0. 15,  and  a  typical  a  is  about  0.5,  which  would  yield 
a  B  of  less  than  10  percent. 

In  conclusion,  the  nonisotropic  flux  distribution  appears  to  be  less  penetrating 
than  an  isotropic  flux  for  a  uniform  shield  thickness.  When  only  a  small  area  is 
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thinly  shielded,  a  nonisotropic  flux  could  create  a  higher  penetrating  dose  than  an 
isotropic  distribution.  Each  shielding  code  contains  parameters  to  describe  the  non¬ 
isotropic  flux,  but  in  some  instances  experimental  data  are  required  to  establish  the 
relationship. 


157 


8.0  PROGRAM  INTEGRATION  AND  UTILIZATION 


The  program  has  been  designed  so  that  the  various  codes  can  be  utilized 
independently.  This  has  advantages  when  studying  the  variations  in  the  environment 
when  no  specific  vehicle  is  given,  or  when  studying  the  effects  of  varying  shielding 
configurations  that  are  independent  of  the  detailed  variations  of  the  environment  or 
of  mission  trajectories.  The  use  of  these  codes  in  this  manner  naturally  requires  the 
analyst  to  exercise  judgment.  This  is  an  essential  element  in  increasing  the  flexi¬ 
bility  of  the  program. 

Once  a  particular  vehicle  has  been  chosen,  it  is  useful  to  automate  the  whole 
trajectory  dose-rate  study.  The  sectoring  code  and  environmental  codes  provide  the 
necessary  inputs  for  the  shielding  code.  The  shielding  code  generates  the  dose 

sach  spectral  type  and  particle.  All 
sec  .  This  dose  conversion  factor 
j  body  point  for  each  particle  type 
are  submitted  as  input  to  the  trajectory  code.  The  trajectory  code  then  automatically 
calculates  the  flux  of  each  spectral  type,  converts  to  dose  rate  and  integrates  the 
flux  and  dose  rate  along  each  of  a  family  of  specified  trajectories.  A  description  of 
how  this  is  accomplished,  using  a  variable  step  integration  routine,  follows. 

The  environment  has  been  placed  on  B-L  grids  in  such  a  way  that  through 
iogarithmic  interpolation  of  the  flux,  and  linear  interpolation  in  the  C  table,  the 
differential  energy  flux  is  obtained  in  the  form: 


delivered  to  each  body  point  per  unit  flux  of 
dose  rates  are  in  rads/hr  from  one  particle  cm 
labeled  k.(j)  for  the  i"1  spectral  type  and  the 


8 

f(B,  L,  E)  *  F(B,  L)  H  C.(B,  L)  f.(E) 
i«l 


The  trajectory  code  gives  the  trajectory  at  any  required  number  of  points  t  ,  prin¬ 
cipally  as  B  (t  ),  L  (t  ).  The  flux  at  t  is  F  (B  ,  L  );  the  differential  fluxnat  t 
.  .  .nnnn  n  nn  n 

is  given  by: 


8 


f<VLn'E)'F(Bn'ln)HC!(8„'Ln)fi(E) 

n=  1 


The  shielding  code  will  predicta  perrunit-particle  dose  conversion,  designated 
by  k.(i),  at  the  j™  body  point  due  to  the  i  basic  spectrum.  Then  the  dose  rate  to 
the  'j™  point  at  time  t  is  given  by 
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d  (t  )=F(B  ,L  )  7“C  (B  ,L  )  k  (j) 

I  n  n  n  * -  inn  i 

i-1 


The  cumulative  flux  to  t  encountered  in  the  mission  is  given  by  the  integral 


F(tn)  -  ^  FtB(t),  L(t)]  dt 


The  dose  encountered  by  the  j^1  point  up  to  t  is  given  by  the  integral 


i  i 

■w-S  diw 

o 

t 

n 

-  ^  FIB(t), L(t)]  YZ C.  CB(t),L(t)3  k.(j)  dt 


o  i 

t 


n 


■ZIsw  ^ 


FCB(t),  L(t)]C.  [B(t),  L(t)3  dt 


The  contribution  to  the  dose  integral  in  a  small  interval  of  time  is  given  by 

8  rn 

AD .frn)-?~ k.(j)  l  Ftt(t),LW)  C,  CB(t),l(t»  dt 

i-1  t  . 

n- 1 

Since  the  coefficient  C.  is  slowly  varying  in  the  interval  we  replace  it  by  its  mean 
valuet 

C(B  L)+C(B  l  ) 

C  (*  \  m  1  n  n _ 1  n  I#  n  I 

iV  2 


but  normalized  to  a  unit  sum; 


8 


y~  c.(t )- 1 

L —  r  n 
i-1 


160 


Then  AD  takes  the  form: 


f0 

AD.(t  )-  <Tk.(i)  C  (t  )  \  FCB(t),L(t)]  dt 

I  n  -  i  i  n  J 


i*  1 


n-1 


,<i) 


(t  )  AF(t  ) 
n  n 


which  involves  the  integral  AF(t),  already  evaluated  in  the  cumulative  flux  integral. 
Effectively  an  average  differential  flux  has  been  defined  for  the  interval  as: 


f(t  ,E) 

n 


**1 

Wi 


Ec,t  )  f. (E) 

L -  in  i 


With  this  assumption  the  integrals  all  reduce  to  an  integral  of  the  flux.  Since  this 
flux  F(B,  L)  is  a  very  strong  function  of  B  and  L  in  some  regions  of  space  and  is  either 
zero  or  a  slow  function  of  B  and  L  in  the  remainder,  a  uniform  time  step  is  not  justi¬ 
fied.  Therefore  we  adopt  a  variable  step  integration  routine  wherein  the  step  is 
controlled  by  the  integrand  —  the  flux. 

Computation  time  will  generally  increase  when  a  variable-step  routine  is 
employed.  To  guard  against  an  unnecessary  increase  in  run  time  due  to  poor  definition 
of  the  flux  map  at  its  edges  or  excessive  use  of  the  C-3  Code,  the  operator  must  be 
provided  with  controls.  A  minimum  time  step  will  eliminate  unnecessary  interval  cut¬ 
backs  in  the  edges  of  the  belt  where  the  flux  is  low  but  rapidly  varying  with  position. 
An  error  criterion  for  the  integral  will  allow  the  optimum  time  step  to  be  chosen  by 
the  routine.  When  the  step  is  such  that  the  error  is  less  than  five  times  the  error 
criterion,  an  8-fold  interpolation  in  B  and  L  will  generally  give  an  acceptable  error. 
Utilizing  the  interpolation  routine  in  this  manner  avoids  excessive  use  of  Mcllwain's 
code.  If  the  step  increases  too  much,  it  is  possible  to  miss  a  significant  flux  at  low 
altitudes  in  the  anomaly.  Controlling  the  time  step  with  a  maximum  allowable  time 
step  will  guard  against  this  possibility.  It  is  recommended  that  a  maximum  time  step 
of  480  sec  be  employed.  Also,  the  variable-step  routine  can  be  eliminated  by 
equating  the  minimum  time  interval  to  the  maximum  time  interval.  Such  a  run  could 
then  be  employed  to  study  the  effectiveness  of  the  variable-step  routine. 

The  following  are  stored  on  magnetic  tape  in  a  format  suitable  for  automatic 
plotting: 
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1)  trajectory  parameters  (t,  B,  L,  R,  X) 

-2  -1 

2)  flux  of  each  spectral  type  In  particle  cm  sec 

_2 

3)  accumulated  flux  of  each  spectral  type  in  particle  cm 

-2  -1 

4)  total  flux  in  particle  cm  sec 

_2 

5)  accumulated  flux  in  particle  cm 

This  tape  provides  the  basic  information  needed  by  the  environment  plotting  code  which 
instructs  a  plotter  to  automatically  display  the  desired  information,  such  as  the  time 
dependence  of  the  dose  rate. 
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9.0  CONCLUSIONS  AND  RECOMMENDATIONS 


9.1  CONCLUSIONS 

It  is  felt  that  the  objectives  of  this  study  have  been  achieved  and  that  a  second- 
generation  computer  program  has  been  developed  commensurate  with  the  experimental 
data  and  theory.  The  program  has  growth  capability  in  directions  presently  indicated 
by  experimental  and  theoretical  studies  and  flexibility  in  its  application  to  conceptual 
and  detailed  design  of  space  systems. 

9. 1. 1  Spacecraft  Trajectories 

Two  methods  have  been  developed  to  treat  the  spacecraft  trajectories.  One 
method  called  approximate  trajectories  computes  the  parameters  r,  <p,  0,  t  from  the 
input  information:  apogee,  altitude,  geographic  location  of  initial  perigee,  and  the 
angle  of  inclination.  A  second  method  will  accept  point-by-point  position  of  the 
spacecraft  from  other  programs,  A  code  has  been  incorporated  to  convert  the  position 
in  geographic  coordinates  to  magnetic  coordinates.  Computer  machine-time-saving 
techniques  have  been  included  to  achieve  a  compromise  between  program  accuracy, 
uncertainties  in  environmental  data,  and  program  cost. 

9. 1. 2  Radiation  Environment 

The  radiation  environment  has  been  arbitrarily  divided  into  trapped  and  untrapped 
radiation  and  the  mathematical  representation  of  these  radiation  environments  has  been 
developed.  This  treatment  includes  energy,  angular,  and  time-dependent  relationships 
for  each  radiation  type.  It  is  believed  that  these  mathematical  representations  will 
handle  any  new  experimental  environmental  data. 

9. 1.3  Vehicle  Analysis 

The  analysis  of  the  effective  shielding  of  the  spacecraft  is  accomplished  by  the 
techniques  called  "sector  analysis. "  In  this  technique  the  dose  point  is  the  origin  of  a 
spherical  coordinate  system  and  the  slant  thicknesses  of  different  materials  for  each 
solid  angle  sector  is  determined.  The  self-shielding  of  the  human  body  is  included  to 
determine  the  dose,  dose  rate,  and  LET  spectrum  of  various  points  in  an  astronaut. 

9.1.4  Electron  Shielding 

A  machine  program  has  been  successfully  developed  to  compute  the  LET  and  dose 
from  penetrating  electrons  inside  a  space  vehicle.  The  program  is  based  on  an  empirical 
relationship  derived  from  Monte  Carlo  data.  These  results  have  been  compared  with 
experimental  and  theoretical  results  of  other  investigators.  The  accuracy  of  the 
empirical  fit  is  within  15  percent  for  normal  incidence,  but  errors  as  great  as  a  factor 
of  two  can  occur  for  isotropic  incidence  when  the  shield  thickness  is  close  to  the  range. 


163 


The  bremsstrahlung  code  has  been  modified  to  describe  volume-distributed 
bremsstrahlung  and  multilayered  shields.  The  significant  result  is  that  the  attenuation 
is  not  changed,  but  the  production  in  each  layer  is  proportional  to  the  atomic  number 
of  the  shield. 

9. 1. 5  Proton  Shielding 

Two  computer  programs  have  been  developed  to  give  a  rapid  and  accurate 
evaluation  of  the  absorbed  dose  resulting  from  the  proton  fluxes  in  space.  The  primary 
proton  program  can  be  used  for  evaluating  the  radiation  shielding  offered  by  the 
majority  of  present  and  planned  space  vehicles.  The  secondary  radiations  are  of 
importance  only  for  fluxes  of  very  soft  protons  or  for  vehicles  with  minimum  shield 
thicknesses  greater  than  about  30  gram/cmz  of  aluminum.  As  thicknesses  greater 
than  60  gram/cm^  are  considered,  the  neutron  dose  becomes  more  important  than  the 
incident  proton  dose.  The  secondary  proton  code  provides  an  upper  limit  for  the 
neutron  dose  contribution.  If  a  more  accurate  evaluation  of  the  neutron  dose  becomes 
necessary,  the  source  terms  of  the  secondary  proton  code  can  provide  a  starting  point 
for  a  more  complete  analysis. 

9. 1. 6  Effect  of  Angular  Distribution  of  Incident  Radiation 

A  study  was  made  to  indicate  the  difference  in  the  internal  flux  if  the  incident 
flux  was  nonisotropic.  All  cases  studied  assumed  that  the  specified  incident  flux  was 
symmetric  about  an  axis.  Three  types  of  particle  deflections  were  studied:  (1)  trans¬ 
mission  where  an  incident  particle  is  degraded  in  energy  but  retains  its  original  direc¬ 
tion,  (2)  transmission  where  the  incident  particle  is  degraded  and  the  distribution  in 
final  direction  is  given  by  f(0),  and  (3)  transmission  of  secondary  particles  with  a 
definite  directional  distribution. 

It  is  concluded  that  the  nonisotropic  flux  distribution  appears  to  be  less 
penetrating  than  an  isotropic  flux  for  a  uniform  cylindrical  shield.  When  only  a  small 
area  of  the  cylinder  is  thinly  shielded,  a  nonisotropic  flux  could  create  a  higher  dose 
than  an  isotropic  distribution.  Each  shielding  code  contains  parameters  to  describe 
the  nonisotropic  flux,  but  in  some  instances  experimental  data  are  required  to 
establish  the  relationship. 

9.2  RECOMMENDATIONS 

Considerable  effort  has  been  made  to  devise  a  second-generation  computer  pro¬ 
gram  with  maximum  growth  capability  and  with  maximum  flexibility;  nevertheless, 
there  are  inevitable  shortcomings  because  of  study  program  limitations,  of  the  rapid 
pace  of  space  research,  and  of  continued  research  in  particle  interactions.  The 
following  specific  areas  of  investigation  are  recommended  to  keep  this  computer  pro¬ 
gram  current: 
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1)  The  time  dependence  of  natural  and  artificially  injected  geomagnetically 
trapped  radiation; 

2)  The  acceptance  cone  method  of  determining  the  magnetic  shielding  of  solar 
particle  radiation; 

3)  Automatic  sectoring  of  vehicle  and  man; 

4)  Extension  of  computer  programs  to  include  primary,  secondary,  and  LET 
spectrum  for  alpha  and  heavier-Z  particles; 

5)  Experimental  verification  of  particle  interaction  theory,  multislab  approxi¬ 
mations,  and  sectoring  approximations; 

6)  The  updating  of  flux  and  spectrum  maps  as  more  current  information 
becomes  available;  and 

7)  Secondary  and  tertiary  components  of  the  LET  spectrum. 
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APPENDIX 


SURVEY  OF  MULTIPLE-SCATTERING  THEORIES 


I.  INTRODUCTION 

The  difficulties  inherent  in  the  problem  of  multiple  scattering  of  fast  charged 
particles  have  given  rise  to  a  diversity  of  theories  and  techniques,  both  analytical 
and  numerical.  Most  of  these  techniques  involve  some  form  of  approximation  and  are 
usually  limited  in  their  domain  of  application.  Some  of  the  theories  predict  only  the 
transmission  factor  or  the  angular  distribution,  while  others  may  include  both,  but 
neglect  various  correction  factors.  Many  "exact"  theories  often  cannot  be  reduced  to 
a  form  convenient  for  numerical  evaluation  or  require  so  much  computation  that  they 
are  unsuited  for  practical  application.  In  view  of  these  difficulties,  the  more  useful 
multiple-scattering  theories  will  be  briefly  reviewed  here  in  order  to  present  the  rela¬ 
tive  merits  of  each  in  terms  of  correlations  with  available  experimental  results  and  the 
case  with  which  computations  may  be  performed. 

II.  ANALYTICAL  METHODS 

The  various  analytic  formulations  of  multiple  scattering  are  based  either  on  the 
random  flight  concept  or  on  the  process  of  diffusion.  The  random-flight  treatment  is 
the  most  rigorous  since  it  exactly  duplicates  the  physical  process  involved  in  multiple 
scattering. 

A.  Random  Flight 

Formal  solutions  of  the  random  flight  problem  have  been  obtained  by  Grosjean 
(Reference  99),  Wigner  (Reference  100),  Guth  and  Inonu  (Reference  101),  Wang  and 
Guth  (Reference  102)  and  Breitenberger  (Reference  103).  The  solutions  of  Grosjean 
and  Wigner  give  the  exact  angular  and  spatial  distributions  but  apply  to  an  infinite 
medium  and  also  depend  on  the  actual  number  of  scatterings  suffered  by  the  particles 
involved.  Breitenberger  has  obtained  the  angular  distribution  for  the  case  of  a  beam  of 
particles  normally  incident  on  a  plane  parallel  slab  of  arbitrary  thickness,  with  pro¬ 
visions  for  Including  energy  loss  within  the  medium.  However,  his  results  are  in  the 
form  of  coupled  recurrence  relations  for  the  scattering  function,  which  are  not  directly 
suited  for  simple  computation.  Explicit  forms  for  the  scattering  function  are  obtained 
by  Breitenberger  by  successive  forward  approximations,  and  the  process  of  extending 
the  results  to  arbitrarily  large  angles  is  outlined.  However,  the  numerical  computa¬ 
tion  Involved  even  at  angles  as  low  as  35  degrees  is  excessive;  and,  in  this  range,  the 
results  are  identical  with  those  obtained  through  the  stmpler  diffusion  process. 

Despite  the  lack  of  useful  solutions,  the  random  flight  technique  is  of  interest 
because  it  indicates  a  fallacy,  which  is  usually  not  mentioned,  in  the  alternative 
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diffusion  process.  According  to  Brettenberger  (Reference  103),  the  diffusion  process 
is  based  on  the  assumption  that  the  directional  distribution  of  the  particles,  after  a 
particular  collision,  does  not  depend  on  the  direction  in  which  other  particles  have 
already  left  the  target  in  earlier  collisions.  But  this  condition  of  "stochastic"  inde¬ 
pendence  of  collisions  holds  only  for  small  scattering  angles. 

B.  Diffusion  Process 

The  majority  of  multiple-scattering  theories  are  based  on  the  assumption  that 
the  scattering  process  is  adequately  described  by  ordinary  diffusion.  This  involves  the 
solution  of  the  integro -differential  transport  equation  with  the  desired  boundary  condi¬ 
tions.  Goudschmidt  and  Sauno'erson  (Reference  53)  obtained  an  expression  for  the 
angular  distribution  by  using  a  Legendre  series  expansion.  Aside  from  the  objection 
stated  by  Breitenberger,  their  result  is  considered  valid  for  all  scattering  angles  and 
can  be  used  with  any  appropriate  (convergent)  single -scattering  cross  section.  The 
effects  of  energy  loss  are  included  in  the  theory,  but  spatial  deflections  resulting  from 
multiple  scattering  are  neglected.  A  difficulty  associated  with  their  result  is  the 
extensive  computation  required  In  evaluating  the  slowly  convergent  Legendre  series. 
Spencer  (Reference  104)  has  derived  simple  recurrence  relations  for  the  series  coeffi¬ 
cients,  but  a  relatively  large  number  of  terms  must  be  summed  for  accurate  results. 

At  large  angles  050  degrees),  Berger  (Reference  54)  has  shown  that  nearly  sixty  terms 
are  required  to  provide  convergence  to  six  significant  figures.  A  further  problem  with 
the  Goudschmidt-Saunderson  theory  is  that  It  is  dependent  on  the  actual  path  traversed 
by  the  particles  rather  than  the  thickness  of  the  target. 

Lewis  (Reference  105)  has  applied  the  Goudschmidt-Saunderson  theory  to  the 
case  of  an  infinite  medium,  taking  the  energy  loss  into  account  by  regarding  it  only 
as  a  function  of  the  residual  range  of  the  scattered  particles.  Some  attempts  have  also 
been  made  to  adapt  the  Goudschmidt-Saunderson  theory  to  the  case  of  bounded  media. 
Side!  et  al.  (Reference  106)  investigated  the  multiple  scattering  of  Compton  recoil 
electrons  produced  by  Co^®  gamma  rays  and  compared  their  results  with  the 
Goudschmidt-Saunderson  theory  by  assuming  the  relation  x  *  t/cos9  between  the  path 
length  of  the  scattered  particle  and  the  target  thickness.  Theoretical  values  were  ob¬ 
tained  by  summing  the  first  eight  terms  in  the  Legendre  series  expansion  of  the  angular 
distribution  function.  Their  results  indicated  reasonably  good  qualitative  agreement 
for  their  targets  (~  20  mg/cm2)  at  most  angles  out  to  140  degrees  but  much  poorer 
agreement  for  thick  targets  (~  200  mg/cm*). 

C.  Small -Angle  Approximations 

Because  of  the  above-mentioned  difficulties  associated  with  general  solutions  of 
the  multiple -scattering  problem,  various  approximate  results  have  been  obtained.  The 
approximations  usually  involve  thin  targets,  small  scattering  angles,  and  no  energy 
loss.  A  complete  summary  of  the  many  mathematical  techniques  developed  has  been 
compiled  by  Scott  (Reference  107).  One  of  the  most  useful  theories  is  that  of  Moiiere 
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(References  52  and  108).  The  Mol  fere  theory  is  basically  just  the  small -angle  limit  of 
the  Goudschmidt-Saunderson  theory,  but  does  not  include  the  effects  of  energy  loss. 
Hanson  et  al.  (Reference  58)  have  shown  that  Mol  fere's  theory  agrees  reasonably 
well  with  their  results  from  thin— foil  scattering  experiments,  and  the  relativistic  cor¬ 
rections  introduced  by  NIgam  et  al.  (Reference  109)  have  improved  the  agreement. 
The  results  indicated  that  the  Mol  iere  scattering  distribution  is  valid  out  to  angles  as 
large  as  20  to  30  degrees.  Bethe  (Reference  52)  has  extended  the  Mol  iere  theory  to 
larger  angles  by  means  of  a  correction  factor.  Although  Breltenberger  has  stated  that 
small -angle  theories  cannot  be  legitimately  extended  to  larger  angles  due  to  the  nature 
of  the  initial  small -angle  approximation,  the  recent  calculations  performed  by  Berger 
(Reference  54)  with  the  Goudschmidt-Saunderson  theory  show  excellent  agreement  with 
Bethe's  large-angle  extension. 

Other  small -angle  scattering  theories  have  been  developed,  notably  by 
Wentzel  (Reference  110)  and  Snider  and  Scott  (Reference  111).  However,  most  of 
these  are  equivalent  to  the  Mol  iere  theory  or  involve  slight  modifications,  usually  in 
the  form  of  various  correction  factors.  Fano  (Reference  112),  for  example,  has  in¬ 
cluded  Inelastic  electron  scattering  by  adding  a  term  to  the  scattering  function.  At 
large  angles,  the  correction  Is  equivalent  to  replacing  Z^  by  Z(Z  +  1).  Cooper  and 
Rainwater  (Reference  113)  and  Ter-Mllcayel lan  (Reference  114)  have  performed  exten¬ 
sive  calculations  of  nuclear  size  effects.  They  Involve  corrections  to  the  screening 
function  in  terms  of  a  nuclear  form  factor  dependent  on  the  nuclear  size.  However, 
according  to  Scott  (Reference  107)  the  form  factor  differs  appreciably  from  unity  only 
for  angles  approximately  10^  times  as  large  as  those  for  which  screening  is  important. 
Thus,  nuclear  size  effects  are  apparently  negligible,  at  least  in  the  small -angle 
approximation. 

Ill  NUMERICAL  METHODS 

The  various  numerical  techniques  which  have  been  developed  for  multiple 
scattering  are  all  based  on  the  random  flight  concept.  The  basic  approach  involves 
grouping  of  the  enormous  number  of  collisions  into  a  sufficiently  small  number  of  repre¬ 
sentative  collisions,  then  sampling  the  "condensed"  random  flights  by  a  Monte  Carlo 
process  according  to  an  appropriate  multiple-scattering  theory.  The  techniques  used 
by  Side!  etal.  (Reference  115),  Leiss  etal.  (Reference  50),  Perkins  (Reference  49), 
Mar  (Reference  44),  and  Berger  (References  54  and  116)  are  all  variations  of  this 
basic  idea. 

The  primary  distinguishing  parameters  in  the  various  approaches  are  the  path 
length  between  each  "condensed"  collision,  the  energy  loss,  angular  deflection,  and 
spatial  displacement.  In  terms  of  these  parameters,  the  most  complete  Monte  Carlo 
model  has  been  developed  by  Berger  (Reference  54).  The  models  of  the  other  refer¬ 
ences  mentioned  above  are  for  the  most  part  restricted  to  particular  choices  of  the 
parameters  with  some  variation  allowed  in  a  single  parameter,  usually  path  length. 
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Berger  used  both  the  Mollere  and  Goudschmidt-Saunderson  scattering  functions  and 
determined  the  path  length  according  to  the  Ionization  and  radiation  energy  losses  by 
requiring  that  the  particle  energy  be  decreased  by  a  certain  fraction  after  each  step. 
Calculation  of  the  transverse  and  longitudinal  displacements  were  also  Included.  Mar, 
Sldel  et  al.,  and  Perkins  all  used  the  small-angle  Mol  fere  theory  but  differed  In  their 
selections  of  path  length.  Perkins  determined  path  length  by  choosing  d/E  to  be  some 
specified  value,  usually  0.02.  Mar  also  used  this  method  along  with  that  of  constant 
path  lengths.  Side!  et  al.  calculated  path  lengths  from  the  Motfere  theory  by  apply¬ 
ing  the  condition  XVB~=  40  degrees.  However,  Sidet  etal.  also  imposed  several 
restricting  approximations,  notably  that  the  final  angular  distribution  Is  described  by 
cos^9  and  that  all  electrons  at  a  particular  depth  In  the  target  have  the  same  energy. 
Letss  etal.  used  a  scattering  distribution  given  by  Rossi  (Reference  117),  equivalent  to 
a  first-order  Gaussian  approximation,  and  a  constant  path  length  of  0. 25  cm. 

Transmission  factors  were  calculated  by  most  of  the  above-mentioned  authors, 
but  only  the  Monte  Carlo  models  developed  by  Mar  and  Berger  were  equipped  to  tabu¬ 
late  angular  distributions  of  transmitted  and  backscattered  electrons.  The  transmission 
and  backscattering  results  for  different  materials  and  Incident  angles  all  agreed  within 
a  few  percent,  with  the  exception  of  the  results  obtained  by  Side!  et  al.  at  low 
energies.  The  deviations,  which  amounted  to  as  much  as  10  to  20  percent,  were  ap¬ 
parently  due  to  the  previously  mentioned  approximation  used  by  Sldei  et  al.  However, 
the  good  qualitative  agreement  with  the  more  elaborate  Monte  Carlo  models  of  the 
other  authors  indicates  the  rather  wide  range  of  capability  of  various  simplifying 
approximations.  Also,  transmission  calculations  performed  by  Berger  with  several  dif¬ 
ferent  scattering  models  and  cross  sections  showed  agreement  to  within  a  few  percent. 

The  angular  distributions  of  transmitted  electrons  obtained  by  Mar  and  Berger 
agreed  reasonably  well  with  each  other  and  with  the  experimental  results  of  Frank 
(Reference  57).  For  penetration  depths  greater  than  the  "diffusion  length,  "  defined 
by  Frank,  the  results  of  Mar  (see  Reference  118)  and  Frank  showed  that  the  angular 
distribution  could  be  described  by  a  cosine -squared  law  In  agreement  with  the  earlier 
assumption  made  by  Sldel  et  al.  Berger1!  results  were  better  fitted  by  (0. 717  cos6  + 
cos^Q)  (Reference  59),  which  indicated  a  relatively  greater  number  of  large  angle 
scatterings. 

Mar  and  Berger  have  also  performed  extensive  computations  in  other  aspects  of 
multiple  scattering.  Mar  has  obtained  the  energy  spectrum  for  transmitted  electrons 
and  has  used  both  isotropic  and  fission  -spectrum  sources.  Investigations  of  multiple 
targets  of  different  materials  are  also  In  progress.  Berger  has  calculated  both  electron 
and  positron  transmissions,  energy  dissipation  in  bounded  and  unbounded  media, 
slowing-down  spectra  and  path  length  straggling,  the  effects  of  isotropic  and  internal 
sources,  and  the  penetration  of  protons. 
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CONCLUSIONS 


The  brief  survey  of  multiple -scattering  theories  indicates  that  the  diffusion  pro¬ 
cess  using  the  transport  equation  is,  at  present,  the  only  technique  that  yields  prac¬ 
tical  analytical  solutions.  The  two  most  useful  theories  appear  to  be  that  of 
Goudschmidt-Saunderson  for  thick  targets,  and  the  Moliere-Bethe  approximation  for 
small  angles  and  thin  targets.  The  primary  difficulty  associated  with  the  Goudschmidt- 
Saunderson  theory  is  the  extensive  computation  required. 

The  combination  of  a  diffusion -type  multiple-scattering  theory  with  Monte 
Carlo  random -samp  ling  techniques  appears  to  be  the  most  convenient  method  available 
for  general  calculations  of  problems  involving  multiple  scattering.  The  Monte  Carlo 
approach  effectively  simulates  the  actual  physical  process  and  can,  therefore,  be  used 
to  obtain  nearly  any  type  of  information  about  the  physical  problem.  The  extensive 
calculations  which  have  already  been  performed  with  this  method  have  shown  that 
excellent  agreement  can  be  obtained  with  experimental  results.  The  only  limitations 
appear  to  be  the  accuracy  of  the  particular  multiple -scattering  function  used  and  the 
large  amount  of  computation  involved. 
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AEC  ACTIVITIES 

US  Atomic  Energy  Commission  (Headquarters  Library,  Reports 
Section),  Mail  Station  G-017,  Wash,  DC  205l*5 

Sandia  Corporation  (Technical  Library),  P.  0.  Box  969,  Livermore, 
Calif  9U551 
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Chief,  Division  of  Technical  Information  Extension,  US  Atomic 
Enerby  Commission,  Box  62,  Oak  Ridge,  Tenn  37831 

University  of  California  Lawrence  Radiation  Laboratory 
(Technical  Information  Division),  P.  0.  Box  808,  Livermore, 
Calif  91*551 

Director,  Los  Alamos  Scientific  Laboratory  (Helen  Redman, 
Report  Library),  P.  0.  Box  1663,  Los  Alamos,  NM  8755^ 

Brookhaven  National  Laboratory,  Upton,  Long  Island,  NY 

Oak  Ridge  National  Laboratory  (Tech  Library),  Oak  Ridge,  Tenn 
37831 

Argonne  National  Laboratory  (Tech  Library),  9700  S.  Cass  Ave, 
Argonne,  Ill  60M0 


OTHER 

Langley  Research  Center  (NASA),  ATTN:  Assoc.  Director, 

Langley  Station,  Hampton,  Va 

National  Bureau  of  Standards,  'Radiological  Equipment  Section, 
Wash  25,  DC 

Central  Intelligence  Agency,  2U30  E  Street  NW,  Wash,  DC  20505 

National  Aeronautics  and  Space  Administration,  Manned  Spacecraft 
Center,  ATTN:  Chief,  Technical  Info  Div,  Houston  1,  Tex 

0TS,  Deaprtment  of  Commerce,  Wash  25,  DC 

Institute  for  Defense  Analysis,  Room  2B257,  The  Pentagon, 

Wash  25,  DC  THRU:  ARPA 

Massachusetts  Institute  of  Technology,  Lincoln  Laboratory 
(Document  Library),  P.  0.  Box  73,  Lexington,  Mass  02173 

Aerospace  Corporation,  P.  0.  Box  95085,  Los  Angeles  1*5 »  Calif 

The  Dikevood  Corporation,  4805  Menaul  Blvd  NE,  Albuquerque,  NM 

Union  Carbide  Research  Institute,  ATTN:  Dr.  Clay  Zerby, 

P.  0.  Box  278,  Tarrytovn,  NY 

Boeing  Aircraft  Company,  ATTN:  Dr.  Bryan  Mar,  P.  0.  Box  3707, 
Seattle  24,  Wash 

General  Dynamics,  ATTN:  Richard  A.  Miller,  P.  0.  Box  748, 

Port  Worth  1,  Texas 

General  Dynamics  Corp,  Nuclear  Instrumentation  Division  (NARF), 
ATTN:  Dr.  Kelley  R.  Spearman,  Fort  Worth  1,  Texas 

S  &  ID,  ATTN:  E.  C.  Olivares,  Dept  093,  Downey,  Calif 
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ASD,  AT*HI:  Mr.  Loren  Pittman,  Wright-Patterson  AFB,  Ohio  45^33 

US  Naval  School  of  Aviation  Medicine,  ATTN:  Dr.  H.  J.  Schaefer, 
Pensacola,  Fla 

SSD  (SSTA,  Lt  Col  E.  D.  Harney),  AFUPO,  Los  Angeles,  Calif  90045 

McDonnel  Aircraft  Corp,  ATTN:  Tom  Galbraith,  Dept  1*01,  St  Louis,  Mo 

(tones  Research  Center,  ATTN:  Dr.  Rudolf  A.  Hofftoan,  Biosatellite 
Project,  Moffett  Field,  Mountain  View,  Calif 

Hughes  Research  Laboratories,  ATTN:  Dr.  Norman  Bailey,  3011 
Malibu  Canyon  Road,  Malibu,  Calif 

Hq,  NASA,  ATTN:  Mr.  Warren  Keller,  1520  H  Street  NW,  Wash  25,  DC 

AVC0  Corporation,  ATTN:  Mr.  Kay  Hoalst,  P.  0.  Drawer  N,  Admiral 
Station,  Tulsa  15 ,  Okla 

S  &  ID,  North  American  Aviation,  ATTN:  Mr.  Fred  Raymes,  Dept 
695/224,  12214  Lakewood  Blvd,  Downey,  Calif 

Bendix  Corp,  Bendix  Systems  Div,  ATTN:  Mr.  Lynn  Lewis,  3300 
Plymouth  Road,  Ann  Arbor,  Mich 

Iowa  State  University,  ATTN:  Dr.  James  Van  Allen,  Dept  of 
Physics  &  Astronomy ,  Ames,  la 

Hq  AFSC  (SCT,  MaJ  R.  E.  Linkous),  Andrews  AFB,  Wash,  DC  20331 

Bell  Laboratories,  ATTN:  Dr.  Walter  L.  Brown,  Murray  Hill,  NJ 

AF0AR  (RR0SE,  Lt  Col  Howard),  Bldg  T-D,  Wash,  DC  20333 

AFSC  (SCGS,  Col  S.  E.  Lifton),  Andrews  AFB,  Wash,  DC  20331 

North rup  Space  Laboratory,  A TTN:  Mr.  Robert  Fortney,  Nuclear 
Sciences  Laboratory,  1111  E.  Broadway,  Hawthorne,  Calif 

Aerospace  Corp,  ATTN:  Dr.  Stan  Freden,  Space  Physics  Laboratory, 

El  Segundo,  Calif 

Aerospace  Corp,  ATTN:  Bob  Pruett,  Los  Angeles,  Calif 

US ATS AM  (SMBN,  Lt  Col  Ed  Ballinger),  Brooks  AFB  Texas  7823? 

Lockheed  Missiles  &  Space  Co.,  ATTN:  Mr.  Bill  Imhof,  Phys 
Sciences  Lab,  Palo  Alto,  Calif 

SSD  (SSTRS,  MaJ  Don  Evans),  AFUPO,  Los  Angeles,  Calif  90045 

NASA  Manned  Space  Flight  Center,  ATTN:  Mr.  Jerry  Modisette, 

Houston  1,  Texas 

AF08R  (SRPH,  Lt  Col  Ralph  Allen),  Bldg  T-D,  Wash,  DC  20333 

General  Dynamics /Astronautics ,  ATTN:  A.  W.  McReynolds,  Space 
Science  Leb  (598-0),  San  Diego,  Calif 
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DISTRIBUTION  (cont'd) 

Dept  of  Defense  (DDR  4  E,  Lt  Col  Lou  Allen),  Wash  25,  DC 

General  Dynamics/Astronautics,  ATTN:  Edward  L.  Noon  (528-40), 

5001  Kearny  Villa  Road,  San  Diego,  Calif 

Hq  USAF  ( AFRST,  Col  Jack  Bollerud) ,  Wash,  DC  20330 

LASL,  ATTN:  Dr.  Wright  Laugh am,  Los  Alamos,  NM 

General  Electric  Co,  ATTN:  Dr.  R.  T.  Frost,  Mgr,  Space  Physics, 

Rm  M1316,  VF  STC,  P.  0.  Box  8555,  Philadelphia  1,  Pa 

Lovelace  Clinic,  ATTN:  Dr.  William  R.  Lovelace,  Albuquerque,  NM 

UCLA,  ATTN:  Dr.  Jerry  McDonald,  Los  Angeles,  Calif 

Lovelace  Clinic,  ATTN:  Dr.  Clayton  S.  White,  Albuquerque,  NM 

LASL,  ATTN:  Dr.  Payne  Harris,  Los  Alamos,  NM 

University  of  California  Lawrence  Radiation  Laboratory,  ATTN: 

Dr.  Harry  West,  Livermore ,  Calif 

ORNL,  ATTN:  Mr.  Ewerett  Blizzard,  Oak  Ridge,  Tenn 

University  of  California,  ATTN:  Dr.  C.  H.  Sondhaus ,  Donner 
Lab,  Livermore,  Calif 

LL0AR  (Col  Hill),  AFUPO,  Los  Angeles  45,  Calif 

University  of  California,  ATTN:  Dr.  Cornelius  A.  Tobias,  Donner 
Lab,  Livermore,  Calif 

SSD  (SSZB,  Lt  Col  T.  Martin),  AFUPO,  Los  Angeles,  Calif  90045 

Goddard  Space  Flight  Center,  NETS:  Dr.  Wilmot  Hess,  Green  Belt,  Md 

AMD  (AMR,  Col  John  Pickering),  Brooks  AFB,  Tex  78235 

AFCRL  (CRZI,  Mr.  Lou  Katz),  Hanscom  Fid,  Bedford,  Mass  01731 

AFCRL,  ATTN:  Mr.  Charles  Reynolds,  Hanscom  Fid,  Bedford,  Mass  01731 

Hq  NASA,  ATTN:  Art  Reetz,  1520  H  Street,  NW,  Wash  25,  DC 

Lovelace  Clinic,  ATTN:  Dr.  A.  H.  Schwichtenberg,  Albuquerque,  NM 

General  Atomic,  A  Division  of  General  Dynamics  Corp,  P.  0.  Box 
608,  San  Diego  12,  Calif 

(Mr.  Roger  Jurgovan) 

(Dr.  Clyde  Jupiter) 

NASA,  1520  H  Street  NW,  Wash,  DC  20546 
(Dr.  Fred  Casal,  Code  MT) 

(Dr.  Leo  Fox,  Code  KBH) 
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DISTRIBUTION  (cont'd) 

NASA,  Marshall  Space  Flight  Center,  Huntsville,  Ala 
(Dr.  R.  D.  Shelton,  Code  M-RP-N) 

(Dr.  Nat  Edmonson) 

(Mr.  Henry  Stearn) 

(Mr.  Martin  Burrell) 

(Mr.  Gus  Whiton) 

NASA,  Goddard  Space  Flight  Center,  ATTN:  Dr.  George  F.  Pieper, 

Greenbelt,  Md 

Langley  Research  Center,  Langley  Station,  Hampton,  Va 
(Dr.  Trutz  Foelsche) 

(Mr.  Harry  Heyson) 

(Dr.  John  Duberg) 

(Mr.  Raymond  Hansen) 

(Mr.  Wayne  Scott) 

(Mr.  William  Hulten) 

NASA  Manned  Spacecraft  Center,  ATTN:  Mr.  Lou  Vosteen,  Houston,  Tex 

Levis  Research  Center,  Cleveland,  Ohio 
(Mr.  Irving  Karp) 

(Mr.  Robert  Hildebrand) 

(Mr.  R.  V.  Brown) 

Jet  Propulsion  Laboratory,  Pasadena,  Calif 
(Mr.  Dvain  F.  Spencer) 

(Dr.  Robert  MacKin) 

Official  Record  Copy  (Lt  D.  A.  Adams,  WLRB-l) 
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